source("script_02_create_objects_for_graphs.R")
library(foreign)
library(plyr)
library(arm)

lo.ci <- .025
hi.ci <- .975
shading<-.7

##################################################################################################################
##PRESIDENTIAL SELECTION
##################################################################################################################
#OWN GOALS / FIGURE 3 
#define distance to nom, in terms of btw old and --positive means on correct side of old nominee

x.axis <- 1.2
y.axis.text <- .7
label.text <- 1.2
x.axis.text <- 1.2
#DW
nomd$pres.dist.from.old.to.nom.dw <- ifelse( ( nomd$pres.sim.dw > nomd$j5old.sim.dw  ), nomd$nom.proj.sim.dw  - nomd$j5old.sim.dw,
																			 ifelse( (nomd$pres.sim.dw < nomd$j5old.sim.dw ), nomd$j5old.sim.dw - nomd$nom.proj.sim.dw ,NA))
nom.level.stats$pos.taking.own.goal.distance.lo.dw <- tapply(nomd$pres.dist.from.old.to.nom.dw, nomd$nominee, quantile,probs=c(lo.ci))
nom.level.stats$pos.taking.own.goal.distance.med.dw <-  tapply(nomd$pres.dist.from.old.to.nom.dw, nomd$nominee, quantile,probs=c(.5))
nom.level.stats$pos.taking.own.goal.distance.hi.dw <- tapply(nomd$pres.dist.from.old.to.nom.dw, nomd$nominee, quantile,probs=c(hi.ci))
nom.level.stats$pos.taking.own.goal.mean.prob.dw <- tapply(nomd$pres.dist.from.old.to.nom.dw <0, nomd$nominee,mean)#prob of an own goal

pdf("Figure_3A.pdf", height=7, width=5)#MAKE FIGURE 3A

y <- 1:length(nom.level.stats$pos.taking.own.goal.mean.prob.dw)
x <- sort (nom.level.stats$pos.taking.own.goal.distance.med.dw)
x.lo <- nom.level.stats$pos.taking.own.goal.distance.lo.dw[order(nom.level.stats$pos.taking.own.goal.distance.med.dw)]
x.hi <- nom.level.stats$pos.taking.own.goal.distance.hi.dw[order(nom.level.stats$pos.taking.own.goal.distance.med.dw)]
prob.own.goal.sorted <- nom.level.stats$pos.taking.own.goal.mean.prob.dw [order(nom.level.stats$pos.taking.own.goal.distance.med.dw, decreasing=F)]
nominee.sorted <- nominee.title.dw[order(nom.level.stats$pos.taking.own.goal.distance.med.dw)]

par(mfrow=c(1,1),mar=c(4,4.5,0,.5))
plot(x,y , xlim=c(-.8,.75), pch=19, col="black", xlab="", ylab="",axes=F,type="n")
polygon(x=c(0,0,par()$usr[1], par()$usr[1]),
             y=par()$usr[c(3,4,4,3)],  ## vertical limits of plot
             col=gray(shading),   ## desired color
             border=F)         ## no border
points(x,y,pch=19)

segments(x.lo, y, x.hi, y)
axis(1, at = seq(-.8,.8,.2), 
     labels=c("-.8","-.6",  "-.4", "-.2", 0, ".2", ".4", ".6",".8"), mgp=c(2,.5,0), cex.axis=x.axis)
axis(2, at =y , label=nominee.sorted, mgp=c(2,.5,0), cex.axis=y.axis.text, las=2)
abline(v=0, col="gray90")
mtext("Distance from old median to NOMINEE\n(higher= toward president)", 1, line=3,cex=x.axis.text)
#add probability for nominees less than zero
ok <- x<0
text(x.hi[ok]+.06,y[ok], round(prob.own.goal.sorted[ok],2), col="black")
text(-.4, 43, "A) Own goals,\nNOMINATE\n(in terms\nof nominee)", cex=1.5)

dev.off()

#BAILEY
bailey.nomd$pres.dist.from.old.to.nom.bly <- ifelse( ( bailey.nomd$pres.sim.bly > bailey.nomd$j5old.sim.bly  ), bailey.nomd$nom.proj.sim.bly  - bailey.nomd$j5old.sim.bly,
																			 ifelse( (bailey.nomd$pres.sim.bly < bailey.nomd$j5old.sim.bly ), bailey.nomd$j5old.sim.bly - bailey.nomd$nom.proj.sim.bly ,NA))
bailey.nom.level.stats$pos.taking.own.goal.distance.lo.bly <- tapply(bailey.nomd$pres.dist.from.old.to.nom.bly, bailey.nomd$nominee, quantile,probs=c(lo.ci))
bailey.nom.level.stats$pos.taking.own.goal.distance.med.bly <-  tapply(bailey.nomd$pres.dist.from.old.to.nom.bly, bailey.nomd$nominee, quantile,probs=c(.5))
bailey.nom.level.stats$pos.taking.own.goal.distance.hi.bly <- tapply(bailey.nomd$pres.dist.from.old.to.nom.bly, bailey.nomd$nominee, quantile,probs=c(hi.ci))
bailey.nom.level.stats$pos.taking.own.goal.mean.prob.bly <- tapply(bailey.nomd$pres.dist.from.old.to.nom.bly <0, bailey.nomd$nominee,mean)#prob of an own goal

pdf("Figure_3B.pdf", height=7, width=5) #FIGURE 3B

y <- 1:length(bailey.nom.level.stats$pos.taking.own.goal.mean.prob.bly)
x <- sort (bailey.nom.level.stats$pos.taking.own.goal.distance.med.bly)
x.lo <- bailey.nom.level.stats$pos.taking.own.goal.distance.lo.bly[order(bailey.nom.level.stats$pos.taking.own.goal.distance.med.bly)]
x.hi <- bailey.nom.level.stats$pos.taking.own.goal.distance.hi.bly[order(bailey.nom.level.stats$pos.taking.own.goal.distance.med.bly)]
prob.own.goal.sorted <- bailey.nom.level.stats$pos.taking.own.goal.mean.prob.bly [order(bailey.nom.level.stats$pos.taking.own.goal.distance.med.bly, decreasing=F)]
nominee.sorted <- nominee.title.bly[order(bailey.nom.level.stats$pos.taking.own.goal.distance.med.bly)]

par(mfrow=c(1,1),mar=c(4,4.5,0,0))
plot(x,y , xlim=c(-2,2), pch=19, col="black", xlab="", ylab="",axes=F,type="n")
polygon(x=c(0,0,par()$usr[1], par()$usr[1]),
             y=par()$usr[c(3,4,4,3)],  ## vertical limits of plot
             col=gray(shading),   ## desired color
             border=F)         ## no border
points(x,y,pch=19)
segments(x.lo, y, x.hi, y)
axis(1, at = seq(-2,2,.5),labels=c(-2,"-1.5","-1","-.5","0",".5","1","1.5","2") , mgp=c(2,.5,0), cex.axis=x.axis)
axis(2, at =y , label=nominee.sorted, mgp=c(2,.5,0), cex.axis=y.axis.text, las=2)
abline(v=0, col="gray90")
mtext("Distance from old median to NOMINEE\n(higher= toward president)", 1, line=3,cex=x.axis.text)
#add probability for nominees less than zero
ok <- x<0
text(x.hi[ok]+.16,y[ok], round(prob.own.goal.sorted[ok],2), col="black")
text(-.7, 31, "B) Own goals, Bailey\n(in terms\nof nominee)", cex=1.5)

dev.off()


###################################################
#2) "AGGRESIVE MISTAKES" moaving new median/nominee. farther away from sen. than old median
#court.outcome AND pos.taking--DEFINE IN TERMS OF BOTH NEW MEDIAN AND OLD MEDIAN
###################################################
#DW
#1st, in terms of new median (court.outcome)
nomd$agg.mistake.court.outcome.amount.dw <- abs( nomd$sen.median.sim.dw- nomd$j5new.sim.dw) - abs(  nomd$sen.median.sim.dw- nomd$j5old.sim.dw) 
nom.level.stats$agg.mistake.court.outcome.distance.lo.dw <- tapply(nomd$agg.mistake.court.outcome.amount.dw, nomd$nominee, quantile,probs=c(lo.ci))
nom.level.stats$agg.mistake.court.outcome.distance.med.dw <-  tapply(nomd$agg.mistake.court.outcome.amount.dw, nomd$nominee, quantile,probs=c(.5))
nom.level.stats$agg.mistake.court.outcome.distance.hi.dw <- tapply(nomd$agg.mistake.court.outcome.amount.dw, nomd$nominee, quantile,probs=c(hi.ci))
nom.level.stats$agg.mistake.court.outcome.mean.prob.dw <- tapply(nomd$agg.mistake.court.outcome.amount.dw > 0, nomd$nominee,mean)#prob of an own goal

#DW
pdf("Figure_4C.pdf", height=7, width=6)

y <- 1:length(nom.level.stats$agg.mistake.court.outcome.mean.prob.dw )
x <- sort (nom.level.stats$agg.mistake.court.outcome.distance.med.dw)
x.lo <- nom.level.stats$agg.mistake.court.outcome.distance.lo.dw[order(nom.level.stats$agg.mistake.court.outcome.distance.med.dw)]
x.hi <- nom.level.stats$agg.mistake.court.outcome.distance.hi.dw[order(nom.level.stats$agg.mistake.court.outcome.distance.med.dw)]
prob.agg.mistake.court.outcome.sorted <- nom.level.stats$agg.mistake.court.outcome.mean.prob.dw [order(nom.level.stats$agg.mistake.court.outcome.distance.med.dw, decreasing=F)]
nominee.sorted <- nominee.title.dw[order(nom.level.stats$agg.mistake.court.outcome.distance.med.dw)]
confirmed.sorted <- nom.level.stats$confirmed[order(nom.level.stats$agg.mistake.court.outcome.distance.med.dw)]
prob.text <- .8

par(mfrow=c(1,1),mar=c(4,4.5,0,0))
plot(x,y , xlim=c(-.8,.8), pch=19, col="black", xlab="", ylab="",axes=F, type="n")
polygon(x=c(0,0,par()$usr[2], par()$usr[2]),
             y=par()$usr[c(3,4,4,3)],  ## vertical limits of plot
             col=gray(shading),   ## desired color
             border=F)         ## no border
segments(x.lo, y, x.hi, y)
points(x[confirmed.sorted==1], y[confirmed.sorted==1], pch=19)
points(x[confirmed.sorted==0], y[confirmed.sorted==0], pch=21, bg="white")

axis(1, at = seq(-.8,.8,.2), 
     labels=c("-.8","-.6",  "-.4", "-.2", 0, ".2", ".4", ".6",".8"), mgp=c(2,.5,0), cex.axis=x.axis)
axis(2, at =y , label=nominee.sorted, mgp=c(2,.5,0), cex.axis=y.axis.text, las=2)
abline(v=0, col="gray90")
#add probability for nominees greater than zero
ok <- x>0
text(x.hi[ok]+.08,y[ok], round(prob.agg.mistake.court.outcome.sorted[ok],2), col="black", cex=prob.text)
text(.5, 5.1, "C) Aggressive\nmistakes,\nNOMINATE\n(in terms of\nnew median)", cex=1.3)
text(0,-5,expression(paste("| ", italic(s[m])- italic(j)[5]^1 , " | - | ", italic(s[m])- italic(j)[5]^0, " |"   )), cex=1.3, xpd=NA)

dev.off()

#2nd, agg.mistake in terms of nominee (pos.taking)
#DW
nomd$agg.mistake.pos.taking.amount.dw <- abs( nomd$sen.median.sim.dw- nomd$nom.proj.sim.dw) - abs(nomd$sen.median.sim.dw- nomd$j5old.sim.dw) 
nom.level.stats$agg.mistake.pos.taking.distance.lo.dw <- tapply(nomd$agg.mistake.pos.taking.amount.dw, nomd$nominee, quantile,probs=c(lo.ci))
nom.level.stats$agg.mistake.pos.taking.distance.med.dw <-  tapply(nomd$agg.mistake.pos.taking.amount.dw, nomd$nominee, quantile,probs=c(.5))
nom.level.stats$agg.mistake.pos.taking.distance.hi.dw <- tapply(nomd$agg.mistake.pos.taking.amount.dw, nomd$nominee, quantile,probs=c(hi.ci))
nom.level.stats$agg.mistake.pos.taking.mean.prob.dw <- tapply(nomd$agg.mistake.pos.taking.amount.dw > 0, nomd$nominee,mean)#prob of an own goal

pdf("Figure_4A.pdf", height=7, width=6)

y <- 1:length(nom.level.stats$agg.mistake.pos.taking.mean.prob.dw )
x <- sort (nom.level.stats$agg.mistake.pos.taking.distance.med.dw)
x.lo <- nom.level.stats$agg.mistake.pos.taking.distance.lo.dw[order(nom.level.stats$agg.mistake.pos.taking.distance.med.dw)]
x.hi <- nom.level.stats$agg.mistake.pos.taking.distance.hi.dw[order(nom.level.stats$agg.mistake.pos.taking.distance.med.dw)]
prob.agg.mistake.pos.taking.sorted <- nom.level.stats$agg.mistake.pos.taking.mean.prob.dw [order(nom.level.stats$agg.mistake.pos.taking.distance.med.dw, decreasing=F)]
nominee.sorted <- nominee.title.dw[order(nom.level.stats$agg.mistake.pos.taking.distance.med.dw)]
confirmed.sorted <- nom.level.stats$confirmed[order(nom.level.stats$agg.mistake.pos.taking.distance.med.dw)]

par(mfrow=c(1,1),mar=c(4,4.5,0,0))
plot(x,y , xlim=c(-.8,.8), pch=19, col="black", xlab="", ylab="",axes=F, type="n")
polygon(x=c(0,0,par()$usr[2], par()$usr[2]),
             y=par()$usr[c(3,4,4,3)],  ## vertical limits of plot
             col=gray(shading),   ## desired color
             border=F)         ## no border

segments(x.lo, y, x.hi, y)
points(x[confirmed.sorted==1], y[confirmed.sorted==1], pch=19)
points(x[confirmed.sorted==0], y[confirmed.sorted==0], pch=21, bg="white")

axis(1, at = seq(-.8,.8,.2), 
     labels=c("-.8","-.6",  "-.4", "-.2", 0, ".2", ".4", ".6",".8"), mgp=c(2,.5,0), cex.axis=x.axis)
axis(2, at =y , label=nominee.sorted, mgp=c(2,.5,0), cex.axis=y.axis.text, las=2)
abline(v=0, col="gray90")
#mtext("| Distance from Nominee to senate median | -\n| distance from old median to senate median | ", 1, line=2.5)#, adj=1)
#add probability for nominees greater than zero
ok <- x>0
text(x.hi[ok]+.1,y[ok], round(prob.agg.mistake.pos.taking.sorted[ok],2), col="black")
text(.5, 4.4, "A) Aggressive\nmistakes,\nNOMINATE\n(in terms\nof nominee)", cex=1.3)
text(0,-5,expression(paste("| ",  italic(s[m]) - italic(n), " | - | ", italic(s[m])- italic(j)[5]^0 , " |"   )), cex=1.3, xpd=NA)

dev.off()

#in terms of new median (court.outcome)
bailey.nomd$agg.mistake.court.outcome.amount.bly <- abs( bailey.nomd$sen.median.sim.bly- bailey.nomd$j5new.sim.bly) - abs(  bailey.nomd$sen.median.sim.bly- bailey.nomd$j5old.sim.bly) 
bailey.nom.level.stats$agg.mistake.court.outcome.distance.lo.bly <- tapply(bailey.nomd$agg.mistake.court.outcome.amount.bly, bailey.nomd$nominee, quantile,probs=c(lo.ci))
bailey.nom.level.stats$agg.mistake.court.outcome.distance.med.bly <-  tapply(bailey.nomd$agg.mistake.court.outcome.amount.bly, bailey.nomd$nominee, quantile,probs=c(.5))
bailey.nom.level.stats$agg.mistake.court.outcome.distance.hi.bly <- tapply(bailey.nomd$agg.mistake.court.outcome.amount.bly, bailey.nomd$nominee, quantile,probs=c(hi.ci))
bailey.nom.level.stats$agg.mistake.court.outcome.mean.prob.bly <- tapply(bailey.nomd$agg.mistake.court.outcome.amount.bly > 0, bailey.nomd$nominee,mean)#prob of an own goal

pdf("Figure_4D.pdf", height=7, width=6)

y <- 1:length(bailey.nom.level.stats$agg.mistake.court.outcome.mean.prob.bly )
x <- sort (bailey.nom.level.stats$agg.mistake.court.outcome.distance.med.bly)
x.lo <- bailey.nom.level.stats$agg.mistake.court.outcome.distance.lo.bly[order(bailey.nom.level.stats$agg.mistake.court.outcome.distance.med.bly)]
x.hi <- bailey.nom.level.stats$agg.mistake.court.outcome.distance.hi.bly[order(bailey.nom.level.stats$agg.mistake.court.outcome.distance.med.bly)]
prob.agg.mistake.court.outcome.sorted <- bailey.nom.level.stats$agg.mistake.court.outcome.mean.prob.bly [order(bailey.nom.level.stats$agg.mistake.court.outcome.distance.med.bly, decreasing=F)]
nominee.sorted <- nominee.title.bly[order(bailey.nom.level.stats$agg.mistake.court.outcome.distance.med.bly)]
confirmed.sorted <- bailey.nom.level.stats$confirmed[order(bailey.nom.level.stats$agg.mistake.court.outcome.distance.med.bly)]
prob.text <- .8

par(mfrow=c(1,1),mar=c(4,4.5,0,0))
plot(x,y , xlim=c(-1,1.4), pch=19, col="black", xlab="", ylab="",axes=F, type="n")
polygon(x=c(0,0,par()$usr[2], par()$usr[2]),
             y=par()$usr[c(3,4,4,3)],  ## vertical limits of plot
             col=gray(shading),   ## desired color
             border=F)         ## no border
segments(x.lo, y, x.hi, y)
points(x[confirmed.sorted==1], y[confirmed.sorted==1], pch=19)
points(x[confirmed.sorted==0], y[confirmed.sorted==0], pch=21, bg="white")

axis(1, at = seq(-1,1.25,.25),
  	labels=c("-1","-.75","-.5",  "-.25", 0, ".25", ".5", ".75","1","1.25"), 
mgp=c(2,.5,0), cex.axis=x.axis)
axis(2, at =y , label=nominee.sorted, mgp=c(2,.5,0), cex.axis=y.axis.text, las=2)
abline(v=0, col="gray90")
#mtext("| Distance from new median to senate median | -\n| distance from old median to senate median | ", 1, line=2.5)#, adj=1)
#add probability for nominees greater than zero
ok <- x>0
text(x.hi[ok]+.08,y[ok], round(prob.agg.mistake.court.outcome.sorted[ok],2), col="black", cex=prob.text)
text(.8, 3.5, "D) Aggressive\nmistakes, Bailey \n(in terms of\nnew median)", cex=1.3)
text(0,-3,expression(paste("| ", italic(s[m])- italic(j)[5]^1 , " | - | ", italic(s[m])- italic(j)[5]^0, " |"   )), cex=1.3, xpd=NA)

dev.off()


#now in terms of nominee
bailey.nomd$agg.mistake.pos.taking.amount.bly <- abs( bailey.nomd$sen.median.sim.bly- bailey.nomd$nom.proj.sim.bly) - abs(bailey.nomd$sen.median.sim.bly- bailey.nomd$j5old.sim.bly) 
bailey.nom.level.stats$agg.mistake.pos.taking.distance.lo.bly <- tapply(bailey.nomd$agg.mistake.pos.taking.amount.bly, bailey.nomd$nominee, quantile,probs=c(lo.ci))
bailey.nom.level.stats$agg.mistake.pos.taking.distance.med.bly <-  tapply(bailey.nomd$agg.mistake.pos.taking.amount.bly, bailey.nomd$nominee, quantile,probs=c(.5))
bailey.nom.level.stats$agg.mistake.pos.taking.distance.hi.bly <- tapply(bailey.nomd$agg.mistake.pos.taking.amount.bly, bailey.nomd$nominee, quantile,probs=c(hi.ci))
bailey.nom.level.stats$agg.mistake.pos.taking.mean.prob.bly <- tapply(bailey.nomd$agg.mistake.pos.taking.amount.bly > 0, bailey.nomd$nominee,mean)#prob of an own goal

pdf("Figure_4B.pdf", height=7, width=6)

y <- 1:length(bailey.nom.level.stats$agg.mistake.pos.taking.mean.prob.bly )
x <- sort (bailey.nom.level.stats$agg.mistake.pos.taking.distance.med.bly)
x.lo <- bailey.nom.level.stats$agg.mistake.pos.taking.distance.lo.bly[order(bailey.nom.level.stats$agg.mistake.pos.taking.distance.med.bly)]
x.hi <- bailey.nom.level.stats$agg.mistake.pos.taking.distance.hi.bly[order(bailey.nom.level.stats$agg.mistake.pos.taking.distance.med.bly)]
prob.agg.mistake.pos.taking.sorted <- bailey.nom.level.stats$agg.mistake.pos.taking.mean.prob.bly [order(bailey.nom.level.stats$agg.mistake.pos.taking.distance.med.bly, decreasing=F)]
nominee.sorted <- nominee.title.bly[order(bailey.nom.level.stats$agg.mistake.pos.taking.distance.med.bly)]
confirmed.sorted <- bailey.nom.level.stats$confirmed[order(bailey.nom.level.stats$agg.mistake.pos.taking.distance.med.bly)]

par(mfrow=c(1,1),mar=c(4,4.5,0,0))
plot(x,y , xlim=c(-1,1.4), pch=19, col="black", xlab="", ylab="",axes=F, type="n")
polygon(x=c(0,0,par()$usr[2], par()$usr[2]),
             y=par()$usr[c(3,4,4,3)],  ## vertical limits of plot
             col=gray(shading),   ## desired color
             border=F)         ## no border
segments(x.lo, y, x.hi, y)
points(x[confirmed.sorted==1], y[confirmed.sorted==1], pch=19)
points(x[confirmed.sorted==0], y[confirmed.sorted==0], pch=21, bg="white")

axis(1, at = seq(-1,1.25,.25),
  	labels=c("-1","-.75","-.5",  "-.25", 0, ".25", ".5", ".75","1","1.25"), 
mgp=c(2,.5,0), cex.axis=x.axis)
axis(2, at =y , label=nominee.sorted, mgp=c(2,.5,0), cex.axis=y.axis.text, las=2)
abline(v=0, col="gray90")

#mtext("| Distance from Nominee to senate median | -\n| distance from old median to senate median | ", 1, line=2.5)#, adj=1)
#add probability for nominees greater than zero
ok <- x>0
text(x.hi[ok]+.1,y[ok], round(prob.agg.mistake.pos.taking.sorted[ok],2), col="black")
text(.8, 3.5, "B) Aggressive\nmistakes, Bailey\n(in terms\nof nominee)", cex=1.3)
text(0,-3,expression(paste("| ",  italic(s[m]) - italic(n), " | - | ", italic(s[m])- italic(j)[5]^0 , " |"   )), cex=1.3, xpd=NA)

dev.off()




##########################################################################################
#IS THE PRESIDENT EVER MEDIAN LOCKED
##########################################################################################
#now calculate predicted location of nominee, calculating regime in each sim#
#first, get SENATE REGIMES
#make senator cutpoints to divide A/B and C/D
#DW
nomd$j4j5mid.sim.dw <- (nomd$j4old.sim.dw + nomd$j5old.sim.dw) / 2
nomd$j5j6mid.sim.dw <- (nomd$j5old.sim.dw + nomd$j6old.sim.dw) / 2 

nomd <- within(nomd, {
	sen.regime.dw <- ifelse( sen.median.sim.dw  <  j4j5mid.sim.dw, "A", 
								ifelse( (j4j5mid.sim.dw <= sen.median.sim.dw) & (sen.median.sim.dw < j5old.sim.dw), "B",
								ifelse( (j5old.sim.dw <=  sen.median.sim.dw)  & (sen.median.sim.dw <=  j5j6mid.sim.dw), "C",
								ifelse( sen.median.sim.dw > j5j6mid.sim.dw, "D", NA))))
	})	

####NOW SET UP DIRECT TESTS (Including 3rd and 4th robust regimes)
#FIRST, FOR court.outcome AND NEARLY court.outcome, DESCRIBE REGIMES
nomd <- within(nomd, {
	presreg.court.outcome.dw <- ifelse( above.dw==1 &  j.exit.sim.dw > j5old.sim.dw,   "restoring",
									ifelse( ( above.dw==1 &  j.exit.sim.dw <= j5old.sim.dw & (sen.regime.dw == "A" | sen.regime.dw=="B")),  "gridlock",
									ifelse( above.dw== 1 & ( j.exit.sim.dw  <= j5old.sim.dw) & ( sen.regime.dw  == "C"), "smaller.shift",
									ifelse(above.dw==1 & ( j.exit.sim.dw  <= j5old.sim.dw) & ( sen.regime.dw == "D"), "maxshift",
									ifelse( above.dw==0 &  j.exit.sim.dw < j5old.sim.dw,   "restoring",
									ifelse( ( above.dw==0 &  j.exit.sim.dw >= j5old.sim.dw & (sen.regime.dw == "C" | sen.regime.dw=="D")),  "gridlock",
									ifelse( above.dw== 0 & ( j.exit.sim.dw  >= j5old.sim.dw) & ( sen.regime.dw  == "B"), "smaller.shift",
									ifelse(above.dw==0 & ( j.exit.sim.dw >= j5old.sim.dw) & ( sen.regime.dw == "A"), "maxshift",
									NA))))))))
#now get point prediction for nearly court.outcome
 ncourt.outcome.point.pred.dw <- ifelse(	presreg.court.outcome.dw=="restoring", pres.sim.dw,
										ifelse(	presreg.court.outcome.dw=="gridlock", j5old.sim.dw,
										ifelse(	presreg.court.outcome.dw=="smaller.shift" & above.dw==1, apply(  matrix(c(pres.sim.dw, (2*sen.median.sim.dw-j5old.sim.dw)),ncol=2) , 1, min),						
										ifelse(	presreg.court.outcome.dw=="smaller.shift" & above.dw==0, apply(  matrix(c(pres.sim.dw, (2*sen.median.sim.dw-j5old.sim.dw)),ncol=2) , 1, max),
										ifelse(	presreg.court.outcome.dw=="maxshift",pres.sim.dw, NA)))))			
									
# now define whether nominee is flip or shift within "smaller.shift"
# shift is if president is interior" to 2S-j -- can get his ideal point  in nearly court.outcome and mixed game
# flip  is if president is exterior to to 2S-j -- can only get S's reflection point  in nearly court.outcome and mixed game 3
#NOTE: "FLIP" means "smaller shift" nominations (in terms of FIGURE 2). "SHIFT" means "Maximum Shift nomination"		
	presreg.court.outcome.smaller.shift.type.dw <- 
			ifelse(above.dw==1 & 	presreg.court.outcome.dw=="smaller.shift"& pres.sim.dw <= 2*sen.median.sim.dw-j5old.sim.dw, "shift",
			ifelse(above.dw==1 & 	presreg.court.outcome.dw=="smaller.shift"& pres.sim.dw > 2*sen.median.sim.dw-j5old.sim.dw, "flip",
			ifelse(above.dw==0 & 	presreg.court.outcome.dw=="smaller.shift"& pres.sim.dw >= 2*sen.median.sim.dw-j5old.sim.dw, "shift",
			ifelse(above.dw==0 & 	presreg.court.outcome.dw=="smaller.shift"& pres.sim.dw < 2*sen.median.sim.dw-j5old.sim.dw, "flip", NA))))	
																	
		#NOW DO pos.taking REGIMES
		presreg.pos.taking.dw <- ifelse( above.dw==1 & (sen.regime.dw == "A" | sen.regime.dw=="B"), "gridlock",
									  	ifelse( above.dw==1 & (sen.regime.dw == "C" | sen.regime.dw=="D"), "smaller.shift",
										  ifelse( above.dw==0 & (sen.regime.dw == "C" | sen.regime.dw=="D"), "gridlock",
										  ifelse( above.dw==0 & (sen.regime.dw == "A" | sen.regime.dw=="B"), "smaller.shift", NA ))))
 #get point prediction for pos.taking
  pos.taking.point.pred.dw <- ifelse(presreg.pos.taking.dw=="gridlock", j5old.sim.dw,
										ifelse(	presreg.pos.taking.dw=="smaller.shift" & above.dw==1, apply(  matrix(c(pres.sim.dw, (2*sen.median.sim.dw-j5old.sim.dw)),ncol=2) , 1, min),						
										ifelse(	presreg.pos.taking.dw=="smaller.shift" & above.dw==0, apply(  matrix(c(pres.sim.dw, (2*sen.median.sim.dw-j5old.sim.dw)),ncol=2) , 1, max), NA)))																 
	presreg.pos.taking.smaller.shift.type.dw <- 
			ifelse(above.dw==1 & 	presreg.pos.taking.dw=="smaller.shift" & pres.sim.dw <= 2*sen.median.sim.dw-j5old.sim.dw, "shift",
			ifelse(above.dw==1 & 	presreg.pos.taking.dw=="smaller.shift" & pres.sim.dw > 2*sen.median.sim.dw-j5old.sim.dw, "flip",
			ifelse(above.dw==0 & 	presreg.pos.taking.dw=="smaller.shift" & pres.sim.dw >= 2*sen.median.sim.dw-j5old.sim.dw, "shift",
			ifelse(above.dw==0 & 	presreg.pos.taking.dw=="smaller.shift" & pres.sim.dw < 2*sen.median.sim.dw-j5old.sim.dw, "flip", NA))))	
	})

tapply(nomd$sen.regime.dw, nomd$nominee, table)
tapply(c(nomd$presreg.court.outcome.dw,nomd$sen.regime.dw), c(nomd$nominee,nomd$nominee) , table)
tapply(c(nomd$presreg.pos.taking.dw,nomd$sen.regime.dw), c(nomd$nominee,nomd$nominee) , table)

#make matrix with frequencies of all categories by nominee
library(zoo)
court.outcome.regime.dist.dw  <- data.frame(t(do.call(merge, lapply(tapply(nomd$presreg.court.outcome.dw, nomd$nominee, table), function(x) zoo(c(x), rownames(x))))))
court.outcome.regime.dist.dw <- within(court.outcome.regime.dist.dw, {
	nominee <- rownames(court.outcome.regime.dist.dw)
	court.outcome.regime.modal.dw  <- apply(court.outcome.regime.dist.dw ,1,which.max)
	court.outcome.regime.modal.dw <- ifelse(court.outcome.regime.modal.dw==1, "flip", ifelse(court.outcome.regime.modal.dw==2, "gridlock", 
											 ifelse(court.outcome.regime.modal.dw=="3", "maxshift", ifelse(court.outcome.regime.modal.dw==4,"restoring", NA))))
  court.outcome.regime.max.percent.dw <- apply(cbind(smaller.shift, gridlock, maxshift, restoring),1,max, na.rm=T)  / n.sims
})
head(court.outcome.regime.dist.dw)

#can we test flip or shift prediction?
table(court.outcome.regime.dist.dw$smaller.shift/n.sims>.5)#NO-only 1

pos.taking.regime.dist.dw  <- data.frame(t(do.call(merge, lapply(tapply(nomd$presreg.pos.taking.dw, nomd$nominee, table), function(x) zoo(c(x), rownames(x))))))
pos.taking.regime.dist.dw <- within(pos.taking.regime.dist.dw, {
  nominee <- rownames(pos.taking.regime.dist.dw )
	pos.taking.regime.modal.dw  <- apply(pos.taking.regime.dist.dw ,1,which.max)
	pos.taking.regime.modal.dw <- ifelse(pos.taking.regime.modal.dw==1, "smaller.shift", ifelse(pos.taking.regime.modal.dw ==2, "gridlock",NA))
  pos.taking.regime.max.percent.dw <- apply(cbind(smaller.shift, gridlock),1,max, na.rm=T)  / n.sims
})
table(pos.taking.regime.dist.dw$pos.taking.regime.max.percent.dw>.8)

#separately, define "lower left" in terms of pres theory figure paper for testing 3rd robust prediction
#need to account for president
nomd <- within(nomd, {
	lower.left.dw <- ifelse( (above.dw ==1  & (sen.regime.dw == "A" | sen.regime.dw == "B") &  j.exit.sim.dw <= j5old.sim.dw )  | 
		(above.dw ==0  & (sen.regime.dw == "C" | sen.regime.dw == "D") &  j.exit.sim.dw > j5old.sim.dw) ,1,0)
})

prob.lower.left.dw <-   tapply (nomd$lower.left.dw, nomd$nominee,mean, na.omit=F)				

temp.df <- data.frame(nominee=court.outcome.regime.dist.dw$nominee,court.outcome.regime.modal.dw=court.outcome.regime.dist.dw$court.outcome.regime.modal.dw,
		court.outcome.regime.max.percent.dw=court.outcome.regime.dist.dw$court.outcome.regime.max.percent.dw,
	 pos.taking.regime.modal.dw= pos.taking.regime.dist.dw$ pos.taking.regime.modal.dw, pos.taking.regime.max.percent.dw= pos.taking.regime.dist.dw$pos.taking.regime.max.percent.dw,prob.lower.left.dw=prob.lower.left.dw )
	
#join to nom-data stats
nomd <- join(nomd,temp.df )
table(nomd$court.outcome.regime.modal.dw,nomd$pos.taking.regime.modal.dw )
head(nomd)

#BAiley
bailey.nomd$j4j5mid.sim.bly <- (bailey.nomd$j4old.sim.bly + bailey.nomd$j5old.sim.bly) / 2
bailey.nomd$j5j6mid.sim.bly <- (bailey.nomd$j5old.sim.bly + bailey.nomd$j6old.sim.bly) / 2 

bailey.nomd <- within(bailey.nomd, {
	sen.regime.bly <- ifelse( sen.median.sim.bly  <  j4j5mid.sim.bly, "A", 
								ifelse( (j4j5mid.sim.bly <= sen.median.sim.bly) & (sen.median.sim.bly < j5old.sim.bly), "B",
								ifelse( (j5old.sim.bly <=  sen.median.sim.bly)  & (sen.median.sim.bly <=  j5j6mid.sim.bly), "C",
								ifelse( sen.median.sim.bly > j5j6mid.sim.bly, "D", NA))))
	})	

####NOW SET UP DIRECT TESTS (Including 3rd and 4th robust regimes)
#FIRST, FOR court.outcome AND NEARLY court.outcome, DESCRIBE REGIMES
bailey.nomd <- within(bailey.nomd, {
	presreg.court.outcome.bly <- ifelse( above.bly==1 &  j.exit.sim.bly > j5old.sim.bly,   "restoring",
									ifelse( ( above.bly==1 &  j.exit.sim.bly <= j5old.sim.bly & (sen.regime.bly == "A" | sen.regime.bly=="B")),  "gridlock",
									ifelse( above.bly== 1 & ( j.exit.sim.bly  <= j5old.sim.bly) & ( sen.regime.bly  == "C"), "smaller.shift",
									ifelse(above.bly==1 & ( j.exit.sim.bly  <= j5old.sim.bly) & ( sen.regime.bly == "D"), "maxshift",
									ifelse( above.bly==0 &  j.exit.sim.bly < j5old.sim.bly,   "restoring",
									ifelse( ( above.bly==0 &  j.exit.sim.bly >= j5old.sim.bly & (sen.regime.bly == "C" | sen.regime.bly=="D")),  "gridlock",
									ifelse( above.bly== 0 & ( j.exit.sim.bly  >= j5old.sim.bly) & ( sen.regime.bly  == "B"), "smaller.shift",
									ifelse(above.bly==0 & ( j.exit.sim.bly >= j5old.sim.bly) & ( sen.regime.bly == "A"), "maxshift",
									NA))))))))
#now get point prediction for nearly court.outcome
 ncourt.outcome.point.pred.bly <- ifelse(	presreg.court.outcome.bly=="restoring", pres.sim.bly,
										ifelse(	presreg.court.outcome.bly=="gridlock", j5old.sim.bly,
										ifelse(	presreg.court.outcome.bly=="smaller.shift" & above.bly==1, apply(  matrix(c(pres.sim.bly, (2*sen.median.sim.bly-j5old.sim.bly)),ncol=2) , 1, min),						
										ifelse(	presreg.court.outcome.bly=="smaller.shift" & above.bly==0, apply(  matrix(c(pres.sim.bly, (2*sen.median.sim.bly-j5old.sim.bly)),ncol=2) , 1, max),
										ifelse(	presreg.court.outcome.bly=="maxshift",pres.sim.bly, NA)))))			
									
# now define whether nominee is flip or shift within "smaller.shift"
# shift is if president is interior" to 2S-j -- can get his ideal point  in nearly court.outcome and mixed game
# flip  is if president is exterior to to 2S-j -- can only get S's reflection point  in nearly court.outcome and mixed game 3
#there are only about 300 simulations where shift is met						
	presreg.court.outcome.smaller.shift.type.bly <- 
			ifelse(above.bly==1 & 	presreg.court.outcome.bly=="smaller.shift"& pres.sim.bly <= 2*sen.median.sim.bly-j5old.sim.bly, "shift",
			ifelse(above.bly==1 & 	presreg.court.outcome.bly=="smaller.shift"& pres.sim.bly > 2*sen.median.sim.bly-j5old.sim.bly, "flip",
			ifelse(above.bly==0 & 	presreg.court.outcome.bly=="smaller.shift"& pres.sim.bly >= 2*sen.median.sim.bly-j5old.sim.bly, "shift",
			ifelse(above.bly==0 & 	presreg.court.outcome.bly=="smaller.shift"& pres.sim.bly < 2*sen.median.sim.bly-j5old.sim.bly, "flip", NA))))	
																	
		#NOW DO pos.taking REGIMES
		presreg.pos.taking.bly <- ifelse( above.bly==1 & (sen.regime.bly == "A" | sen.regime.bly=="B"), "gridlock",
									  	ifelse( above.bly==1 & (sen.regime.bly == "C" | sen.regime.bly=="D"), "smaller.shift",
										  ifelse( above.bly==0 & (sen.regime.bly == "C" | sen.regime.bly=="D"), "gridlock",
										  ifelse( above.bly==0 & (sen.regime.bly == "A" | sen.regime.bly=="B"), "smaller.shift", NA ))))
 #get point prediction for pos.taking
  pos.taking.point.pred.bly <- ifelse(presreg.pos.taking.bly=="gridlock", j5old.sim.bly,
										ifelse(	presreg.pos.taking.bly=="smaller.shift" & above.bly==1, apply(  matrix(c(pres.sim.bly, (2*sen.median.sim.bly-j5old.sim.bly)),ncol=2) , 1, min),						
										ifelse(	presreg.pos.taking.bly=="smaller.shift" & above.bly==0, apply(  matrix(c(pres.sim.bly, (2*sen.median.sim.bly-j5old.sim.bly)),ncol=2) , 1, max), NA)))																 
	presreg.pos.taking.smaller.shift.type.bly <- 
			ifelse(above.bly==1 & 	presreg.pos.taking.bly=="smaller.shift" & pres.sim.bly <= 2*sen.median.sim.bly-j5old.sim.bly, "shift",
			ifelse(above.bly==1 & 	presreg.pos.taking.bly=="smaller.shift" & pres.sim.bly > 2*sen.median.sim.bly-j5old.sim.bly, "flip",
			ifelse(above.bly==0 & 	presreg.pos.taking.bly=="smaller.shift" & pres.sim.bly >= 2*sen.median.sim.bly-j5old.sim.bly, "shift",
			ifelse(above.bly==0 & 	presreg.pos.taking.bly=="smaller.shift" & pres.sim.bly < 2*sen.median.sim.bly-j5old.sim.bly, "flip", NA))))	
	})

tapply(bailey.nomd$sen.regime.bly, bailey.nomd$nominee, table)
tapply(c(bailey.nomd$presreg.court.outcome.bly,bailey.nomd$sen.regime.bly), c(bailey.nomd$nominee,bailey.nomd$nominee) , table)
tapply(c(bailey.nomd$presreg.pos.taking.bly,bailey.nomd$sen.regime.bly), c(bailey.nomd$nominee,bailey.nomd$nominee) , table)

#make matrix with frequencies of all categories by nominee
library(zoo)
court.outcome.regime.dist.bly  <- data.frame(t(do.call(merge, lapply(tapply(bailey.nomd$presreg.court.outcome.bly, bailey.nomd$nominee, table), function(x) zoo(c(x), rownames(x))))))
court.outcome.regime.dist.bly <- within(court.outcome.regime.dist.bly, {
	nominee <- rownames(court.outcome.regime.dist.bly)
	court.outcome.regime.modal.bly  <- apply(court.outcome.regime.dist.bly ,1,which.max)
	court.outcome.regime.modal.bly <- ifelse(court.outcome.regime.modal.bly==1, "flip", ifelse(court.outcome.regime.modal.bly==2, "gridlock", 
											 ifelse(court.outcome.regime.modal.bly=="3", "maxshift", ifelse(court.outcome.regime.modal.bly==4,"restoring", NA))))
  court.outcome.regime.max.percent.bly <- apply(cbind(smaller.shift, gridlock, maxshift, restoring),1,max, na.rm=T)  / n.sims
})
head(court.outcome.regime.dist.bly)

#can we test flip or shift prediction?
table(court.outcome.regime.dist.bly$smaller.shift/n.sims>.5)#NO-only 2

pos.taking.regime.dist.bly  <- data.frame(t(do.call(merge, lapply(tapply(bailey.nomd$presreg.pos.taking.bly, bailey.nomd$nominee, table), function(x) zoo(c(x), rownames(x))))))
pos.taking.regime.dist.bly <- within(pos.taking.regime.dist.bly, {
  nominee <- rownames(pos.taking.regime.dist.bly )
	pos.taking.regime.modal.bly  <- apply(pos.taking.regime.dist.bly ,1,which.max)
	pos.taking.regime.modal.bly <- ifelse(pos.taking.regime.modal.bly==1, "smaller.shift", ifelse(pos.taking.regime.modal.bly ==2, "gridlock",NA))
  pos.taking.regime.max.percent.bly <- apply(cbind(smaller.shift, gridlock),1,max, na.rm=T)  / n.sims
})
table(pos.taking.regime.dist.bly$pos.taking.regime.max.percent.bly>.8)

#separately, define "lower left" in terms of pres theory figure paper for testing 3rd robust prediction
#need to account for president
bailey.nomd <- within(bailey.nomd, {
	lower.left.bly <- ifelse( (above.bly ==1  & (sen.regime.bly == "A" | sen.regime.bly == "B") &  j.exit.sim.bly <= j5old.sim.bly )  | 
		(above.bly ==0  & (sen.regime.bly == "C" | sen.regime.bly == "D") &  j.exit.sim.bly > j5old.sim.bly) ,1,0)
})

prob.lower.left.bly <-   tapply (bailey.nomd$lower.left.bly, bailey.nomd$nominee,mean, na.omit=F)				

temp.df <- data.frame(nominee=court.outcome.regime.dist.bly$nominee,court.outcome.regime.modal.bly=court.outcome.regime.dist.bly$court.outcome.regime.modal.bly,
		court.outcome.regime.max.percent.bly=court.outcome.regime.dist.bly$court.outcome.regime.max.percent.bly,
	 pos.taking.regime.modal.bly= pos.taking.regime.dist.bly$ pos.taking.regime.modal.bly, pos.taking.regime.max.percent.bly= pos.taking.regime.dist.bly$pos.taking.regime.max.percent.bly,prob.lower.left.bly=prob.lower.left.bly )
	
#join to nom-data stats
bailey.nomd <- join(bailey.nomd,temp.df )
table(bailey.nomd$court.outcome.regime.modal.bly,bailey.nomd$pos.taking.regime.modal.bly )
head(bailey.nomd)

#take nominations where predicted to be in gridlock across all regimes
#compare actual versus j5ol

#NOMINATE
nomd$gridlock.robust <- ifelse(nomd$presreg.court.outcome.dw=="gridlock" & nomd$presreg.pos.taking.dw=="gridlock","gridlock","not gridlock")
gridlock.robust.dw  <- data.frame(t(do.call(merge, lapply(tapply(nomd$gridlock.robust , nomd$nominee, table), function(x) zoo(c(x), rownames(x))))))
gridlock.robust.dw$gridlock <- ifelse(is.na(gridlock.robust.dw$gridlock), 0, gridlock.robust.dw$gridlock)
gridlock.robust.dw$not.gridlock <- ifelse(is.na(gridlock.robust.dw$not.gridlock), 0, gridlock.robust.dw$not.gridlock)

gridlock.robust.dw <- within(gridlock.robust.dw, {
	nominee <- rownames(gridlock.robust.dw)
	gridlock.robust.modal.dw  <- apply(gridlock.robust.dw ,1,which.max)
	gridlock.robust.modal.dw <- ifelse(gridlock.robust.modal.dw==2, "not gridlock", ifelse(gridlock.robust.modal.dw==1, "gridlock", NA))
  gridlock.robust.max.percent.dw <- apply(cbind(gridlock),1,max, na.rm=T)  / n.sims
})
head(gridlock.robust.dw)
nomd <- join(nomd,gridlock.robust.dw)

#Bailey
bailey.nomd$gridlock.robust <- ifelse(bailey.nomd$presreg.court.outcome.bly=="gridlock" & bailey.nomd$presreg.pos.taking.bly=="gridlock","gridlock","not gridlock")
gridlock.robust.bly  <- data.frame(t(do.call(merge, lapply(tapply(bailey.nomd$gridlock.robust , bailey.nomd$nominee, table), function(x) zoo(c(x), rownames(x))))))
gridlock.robust.bly$gridlock <- ifelse(is.na(gridlock.robust.bly$gridlock), 0, gridlock.robust.bly$gridlock)
gridlock.robust.bly$not.gridlock <- ifelse(is.na(gridlock.robust.bly$not.gridlock), 0, gridlock.robust.bly$not.gridlock)

gridlock.robust.bly <- within(gridlock.robust.bly, {
	nominee <- rownames(gridlock.robust.bly)
	gridlock.robust.modal.bly  <- apply(gridlock.robust.bly ,1,which.max)
	gridlock.robust.modal.bly <- ifelse(gridlock.robust.modal.bly==2, "not gridlock", ifelse(gridlock.robust.modal.bly==1, "gridlock", NA))
  gridlock.robust.max.percent.bly <- apply(cbind(gridlock),1,max, na.rm=T)  / n.sims
})
head(gridlock.robust.bly)

bailey.nomd <- join(bailey.nomd,gridlock.robust.bly)

######################################################################################
##MAKE GRAPH for Figure 5
#################################################################################################################################
threshold.hi <-.5

pdf("Figure_5.pdf", height=10, width=12)

x.axis <-1.2
y.axis.text <-1.1
label.text <- 1.2

#DW
ok <- nomd$gridlock.robust.max.percent.dw >= threshold.hi
temp.grid.data.dw <- subset(nomd, ok)#make temp data
#calculate confidence intervals and median estimate
temp.grid.data.dw$nom.v.old <- temp.grid.data.dw$nom.proj.sim.dw-temp.grid.data.dw$j5old.sim.dw 

nom.v.old.lo.dw <- tapply(temp.grid.data.dw$nom.v.old, temp.grid.data.dw$nominee, quantile,probs=c(lo.ci))
nom.v.old.med.dw <-  tapply(temp.grid.data.dw$nom.v.old, temp.grid.data.dw$nominee, quantile,probs=c(.5))
nom.v.old.hi.dw <- tapply(temp.grid.data.dw$nom.v.old, temp.grid.data.dw$nominee, quantile,probs=c(hi.ci))
nom.v.old.mean.prob.dw <- tapply(temp.grid.data.dw$nom.v.old > 0, temp.grid.data.dw$nominee,mean)

dem.pres.temp <- ifelse(tapply(temp.grid.data.dw$pres.sim.dw<0,temp.grid.data.dw$nominee, mean)>0,1,0)
confirmed.temp <-  tapply(temp.grid.data.dw$confirmed > 0, temp.grid.data.dw$nominee,mean)
temp.df <- data.frame(nominee=unique(temp.grid.data.dw$nominee.title.dw),
	nom.v.old.lo.dw,nom.v.old.med.dw,nom.v.old.hi.dw, nom.v.old.mean.prob.dw, dem.pres.temp, confirmed.temp)
temp.df <- temp.df[order(temp.df$dem.pres.temp,nom.v.old.med.dw ),]
	
y <- 1:nrow(temp.df)
x <-temp.df$nom.v.old.med.dw
x.lo <- temp.df$nom.v.old.lo.dw
x.hi <- temp.df$nom.v.old.hi.dw
nominee.sorted <- temp.df$nominee
confirmed.sorted <- temp.df$confirmed.temp
prob.text <- .8

par(mfrow=c(1,2),mar=c(3,6.5,.5,1))
plot(x,y , xlim=c(-.5,.5), pch=19, col="black", xlab="", ylab="",axes=F, type="n")
polygon(x=c(par()$usr[1],par()$usr[1], par()$usr[2],par()$usr[2]), y=c(11.5,par()$usr[4]-.5,par()$usr[4]-.5,11.5), col="gray90")
segments(x.lo, y, x.hi, y)
points(x[confirmed.sorted==1], y[confirmed.sorted==1], pch=19)
points(x[confirmed.sorted==0], y[confirmed.sorted==0], pch=21, bg="white")

axis(1, at = seq(-1,1.4,.2), 
     labels=c("-1","-.8","-.6",  "-.4", "-.2", 0, ".2", ".4", ".6",".8","1","1.2","1.4"), mgp=c(2,.5,0), cex.axis=x.axis)
axis(2, at =y , label=nominee.sorted, mgp=c(2,.5,0), cex.axis=y.axis.text, las=2)
segments(0,par()$usr[4]-.5,0,par()$usr[3], col="gray80")

mtext("Nominee-old median", 1, line=2,cex=label.text)
mtext("A) NOMINATE",3, line=-1, cex=label.text, font=2)

text(-.4, 18, "Democratic\nappointees",cex=label.text, font=3, srt=90)

#Bailey
ok <- bailey.nomd$gridlock.robust.max.percent.bly >= threshold.hi
temp.grid.data.bly <- subset(bailey.nomd, ok)#make temp data
#calculate confidence intervals and median estimate
temp.grid.data.bly$nom.v.old <- temp.grid.data.bly$nom.proj.sim.bly-temp.grid.data.bly$j5old.sim.bly 

nom.v.old.lo.bly <- tapply(temp.grid.data.bly$nom.v.old, temp.grid.data.bly$nominee, quantile,probs=c(lo.ci))
nom.v.old.med.bly <-  tapply(temp.grid.data.bly$nom.v.old, temp.grid.data.bly$nominee, quantile,probs=c(.5))
nom.v.old.hi.bly <- tapply(temp.grid.data.bly$nom.v.old, temp.grid.data.bly$nominee, quantile,probs=c(hi.ci))
nom.v.old.mean.prob.bly <- tapply(temp.grid.data.bly$nom.v.old > 0, temp.grid.data.bly$nominee,mean)

dem.pres.temp <- ifelse( tapply(temp.grid.data.bly$pres.sim.bly<0,temp.grid.data.bly$nominee, mean) > .5,1,0)#need cutoff of .5 to separate
confirmed.temp <-  tapply(temp.grid.data.bly$confirmed > 0, temp.grid.data.bly$nominee,mean)
temp.df <- data.frame(nominee=unique(temp.grid.data.bly$nominee.title.bly),
	nom.v.old.lo.bly,nom.v.old.med.bly,nom.v.old.hi.bly, nom.v.old.mean.prob.bly, dem.pres.temp, confirmed.temp)
temp.df <- temp.df[order(temp.df$dem.pres.temp,nom.v.old.med.bly ),]
	
y <- 1:nrow(temp.df)
x <-temp.df$nom.v.old.med.bly
x.lo <- temp.df$nom.v.old.lo.bly
x.hi <- temp.df$nom.v.old.hi.bly
nominee.sorted <- temp.df$nominee
confirmed.sorted <- temp.df$confirmed.temp
prob.text <- .8

plot(x,y , xlim=c(-1.4,1.4), pch=19, col="black", xlab="", ylab="",axes=F, type="n")
polygon(x=c(par()$usr[1],par()$usr[1], par()$usr[2],par()$usr[2]), y=c(11.5,par()$usr[4]-.4,par()$usr[4]-.4,11.5), col="gray90")
segments(x.lo, y, x.hi, y)
points(x[confirmed.sorted==1], y[confirmed.sorted==1], pch=19)
points(x[confirmed.sorted==0], y[confirmed.sorted==0], pch=21, bg="white")
axis(1, at = seq(-1.4,1.4,.2), 
     labels=c("-1.4","-1.2","-1","-.8","-.6",  "-.4", "-.2", 0, ".2", ".4", ".6",".8","1","1.2","1.4"), mgp=c(2,.5,0), cex.axis=x.axis)
axis(2, at =y , label=nominee.sorted, mgp=c(2,.5,0), cex.axis=y.axis.text, las=2)
segments(0,par()$usr[4]-.5,0,par()$usr[3], col="gray80")
mtext("Nominee-old median", 1, line=2,cex=label.text)
mtext("B) Bailey",3, line=-1, cex=label.text, font=2)
text(-1.2, 13, "Democratic\nappointees",cex=label.text, font=3, srt=90)

dev.off()

##################
#PRESIDENTIAL REGRESSIONS
#DW
#####
#REGRESSIONS OF NOMINEE SELECTION ??????????????????????
#pos.taking
#test is should vary with old median in gridlock, min(p, 2s-j5_old) in flip shift
#use interactions
#define regimes and minimum prediction
#issue: For "Flip", need to figure out if it's p (shift) or flip
nomd <- within (nomd, {
	#court.outcome MODELS
		#CREATE DUMMIES
	court.outcome.gridlock.dw <- ifelse(presreg.court.outcome.dw=="gridlock", 1,0)
	court.outcome.flip.dw <- ifelse(presreg.court.outcome.smaller.shift.type.dw=="flip", 1,0)
	court.outcome.flip.dw <- ifelse( is.na(court.outcome.flip.dw), 0,court.outcome.flip.dw)
	court.outcome.shift.dw <- ifelse(presreg.court.outcome.smaller.shift.type.dw=="shift", 1,0)
	court.outcome.shift.dw <- ifelse( is.na(court.outcome.shift.dw), 0,court.outcome.shift.dw)
	court.outcome.restoring.dw <- ifelse(presreg.court.outcome.dw=="restoring", 1,0)
	court.outcome.maxshift.dw <- ifelse(presreg.court.outcome.dw=="maxshift", 1,0)
	court.outcome.p.combined.dw <-ifelse(	court.outcome.restoring.dw==1 | 	court.outcome.maxshift.dw==1 |	court.outcome.shift.dw== 1 , 1,0)
	court.outcome.p.combined.dw <- ifelse( is.na(court.outcome.p.combined.dw), 0,court.outcome.p.combined.dw)
	#for combining regions
	court.outcome.not.p.combined.dw <- ifelse(	court.outcome.p.combined.dw==0,1,0)
	court.outcome.not.gridlock.dw <- ifelse(	court.outcome.gridlock.dw==0,1,0)
	court.outcome.not.flip.dw <- ifelse(court.outcome.flip.dw==0,1,0)
	#pos.taking MODELS
	pos.taking.gridlock.dw <- ifelse(presreg.pos.taking.dw=="gridlock", 1,0)
	pos.taking.flip.dw <- ifelse(presreg.pos.taking.smaller.shift.type.dw== "flip", 1,0)
	pos.taking.flip.dw <- ifelse( is.na(pos.taking.flip.dw ), 0,pos.taking.flip.dw )
	pos.taking.shift.dw <- ifelse(presreg.pos.taking.smaller.shift.type.dw== "shift", 1,0)
	pos.taking.shift.dw <- ifelse( is.na(pos.taking.shift.dw ), 0,pos.taking.shift.dw )
	pos.taking.not.gridlock.dw <- ifelse(	pos.taking.gridlock.dw==0,1,0)
	pos.taking.not.flip.dw <- ifelse(	pos.taking.flip.dw==0,1,0)
	pos.taking.not.shift.dw <- ifelse(	pos.taking.shift.dw==0,1,0)
	#Senate indifference POINT: this applies across models
	sen.ind.point.dw <- 2*sen.median.sim.dw-j5old.sim.dw	 
})

lo <-.025
hi <- .975

#ncourt.outcome model#
#Model 1: just theory predictions
court.outcome.model.1.dw.intercept <- rep(NA, n.sims)
court.outcome.model.1.dw.gridlock <- rep(NA, n.sims)
court.outcome.model.1.dw.p.predicted <- rep(NA, n.sims)
court.outcome.model.1.dw.flip <- rep(NA, n.sims)
court.outcome.model.1.dw.rsquared <- rep(NA, n.sims)

for (i in 1:n.sims){
	ok <- nomd$simulation==i
  foo <-	lm(nom.proj.sim.dw ~  	court.outcome.gridlock.dw:j5old.sim.dw + court.outcome.p.combined.dw:pres.sim.dw +	court.outcome.flip.dw:sen.ind.point.dw ,	
		data=nomd,subset=ok)
	sim.foo <- sim(foo,n.sims=1)
	court.outcome.model.1.dw.intercept[i] <- coef(sim.foo)[1]
	court.outcome.model.1.dw.gridlock[i] <- coef(sim.foo)[2]
	court.outcome.model.1.dw.p.predicted[i] <- coef(sim.foo)[3]
	court.outcome.model.1.dw.flip[i] <- coef(sim.foo)[4]
	court.outcome.model.1.dw.rsquared[i] <- summary(foo)$r.squared
}		

#MODEL 2: now in non-theory predictions
court.outcome.model.2.dw.intercept <- rep(NA, n.sims)
court.outcome.model.2.dw.not.gridlock <- rep(NA, n.sims)
court.outcome.model.2.dw.not.p.predicted <- rep(NA, n.sims)
court.outcome.model.2.dw.not.flip <- rep(NA, n.sims)
court.outcome.model.2.dw.rsquared <- rep(NA, n.sims)

for (i in 1:n.sims){
	ok <- nomd$simulation==i
  foo <-	lm(nom.proj.sim.dw ~  	
		+ 	court.outcome.not.gridlock.dw:j5old.sim.dw + court.outcome.not.p.combined.dw:pres.sim.dw +	court.outcome.not.flip.dw:sen.ind.point.dw
		, data=nomd,subset=ok)
	sim.foo <- sim(foo,n.sims=1)
	court.outcome.model.2.dw.intercept[i]    <- coef(sim.foo)[1]
	court.outcome.model.2.dw.not.gridlock[i]    <- coef(sim.foo)[2]
	court.outcome.model.2.dw.not.p.predicted[i] <- coef(sim.foo)[3]
	court.outcome.model.2.dw.not.flip[i]        <- coef(sim.foo)[4]	
	court.outcome.model.2.dw.rsquared[i] <- summary(foo)$r.squared
}		
		
############################################
######################
#pos.taking model
######################
######################
#1) just theory predictions
pos.taking.model.1.dw.intercept <- rep(NA, n.sims)
pos.taking.model.1.dw.gridlock <- rep(NA, n.sims)
pos.taking.model.1.dw.shift <- rep(NA, n.sims)
pos.taking.model.1.dw.flip <- rep(NA, n.sims)
pos.taking.model.1.dw.rsquared <- rep(NA, n.sims)
	
for (i in 1:n.sims){
	ok <- nomd$simulation==i
	#pos.taking
 foo <-	lm(nom.proj.sim.dw ~  	pos.taking.gridlock.dw:j5old.sim.dw + pos.taking.shift.dw:pres.sim.dw + pos.taking.flip.dw:sen.ind.point.dw,
		 	data=nomd,subset=ok)
	sim.foo <- sim(foo,n.sims=1)
	pos.taking.model.1.dw.intercept[i] <- coef(sim.foo)[1]
	pos.taking.model.1.dw.gridlock[i] <- coef(sim.foo)[2]
	pos.taking.model.1.dw.shift[i] <- coef(sim.foo)[3]
	pos.taking.model.1.dw.flip[i] <- coef(sim.foo)[4]
	pos.taking.model.1.dw.rsquared[i] <- summary(foo)$r.squared
}		

#2) now non-theory specific predictgions
pos.taking.model.2.dw.intercept <- rep(NA, n.sims)
pos.taking.model.2.dw.not.gridlock <- rep(NA, n.sims)
pos.taking.model.2.dw.not.shift <- rep(NA, n.sims)
pos.taking.model.2.dw.not.flip <- rep(NA, n.sims)
pos.taking.model.2.dw.rsquared <- rep(NA, n.sims)

for (i in 1:n.sims){
	ok <- nomd$simulation==i
	#pos.taking
 foo <-	lm(nom.proj.sim.dw ~  	pos.taking.not.gridlock.dw:j5old.sim.dw + pos.taking.not.shift.dw:pres.sim.dw + pos.taking.not.flip.dw:sen.ind.point.dw, data=nomd,subset=ok)
	sim.foo <- sim(foo,n.sims=1)
	pos.taking.model.2.dw.intercept[i] <- coef(foo)[1]
	pos.taking.model.2.dw.not.gridlock[i] <-  coef(foo)[2]
	pos.taking.model.2.dw.not.shift[i] <- coef(foo)[3]
  pos.taking.model.2.dw.not.flip[i]  <-	 coef(foo)[4]
	pos.taking.model.2.dw.rsquared[i] <- summary(foo)$r.squared
}

#now combined
pos.taking.model.3.dw.intercept <- rep(NA, n.sims)
pos.taking.model.3.dw.gridlock <- rep(NA, n.sims)
pos.taking.model.3.dw.shift <- rep(NA, n.sims)
pos.taking.model.3.dw.flip <- rep(NA, n.sims)
pos.taking.model.3.dw.not.gridlock <- rep(NA, n.sims)
pos.taking.model.3.dw.not.shift <- rep(NA, n.sims)
pos.taking.model.3.dw.not.flip <- rep(NA, n.sims)
pos.taking.model.3.dw.rsquared <- rep(NA, n.sims)

for (i in 1:n.sims){
	ok <- nomd$simulation==i
	#pos.taking
 foo <-	lm(nom.proj.sim.dw ~  	pos.taking.gridlock.dw:j5old.sim.dw     + pos.taking.shift.dw:pres.sim.dw     + pos.taking.flip.dw:sen.ind.point.dw #theory predictions
													+	pos.taking.not.gridlock.dw:j5old.sim.dw + pos.taking.not.shift.dw:pres.sim.dw + pos.taking.not.flip.dw:sen.ind.point.dw, 
													data=nomd,subset=ok)
	sim.foo <- sim(foo,n.sims=1)
	pos.taking.model.3.dw.intercept[i] <- coef(foo)[1]
	pos.taking.model.3.dw.gridlock[i] <- coef(foo)[2]
	pos.taking.model.3.dw.shift[i] <- coef(foo)[3]
	pos.taking.model.3.dw.flip[i] <- coef(foo)[4]
	pos.taking.model.3.dw.not.gridlock[i] <-  coef(foo)[5]
	pos.taking.model.3.dw.not.shift[i] <- coef(foo)[6]
  pos.taking.model.3.dw.not.flip[i]  <-	 coef(foo)[7]
	pos.taking.model.3.dw.rsquared[i] <- summary(foo)$r.squared
}

#BAILEY
bailey.nomd <- within (bailey.nomd, {
	#court.outcome MODELS
		#CREATE DUMMIES
	court.outcome.gridlock.bly <- ifelse(presreg.court.outcome.bly=="gridlock", 1,0)
	court.outcome.flip.bly <- ifelse(presreg.court.outcome.smaller.shift.type.bly=="flip", 1,0)
	court.outcome.flip.bly <- ifelse( is.na(court.outcome.flip.bly), 0,court.outcome.flip.bly)
	court.outcome.shift.bly <- ifelse(presreg.court.outcome.smaller.shift.type.bly=="shift", 1,0)
	court.outcome.shift.bly <- ifelse( is.na(court.outcome.shift.bly), 0,court.outcome.shift.bly)
	court.outcome.restoring.bly <- ifelse(presreg.court.outcome.bly=="restoring", 1,0)
	court.outcome.maxshift.bly <- ifelse(presreg.court.outcome.bly=="maxshift", 1,0)
	court.outcome.p.combined.bly <-ifelse(	court.outcome.restoring.bly==1 | 	court.outcome.maxshift.bly==1 |	court.outcome.shift.bly== 1 , 1,0)
	court.outcome.p.combined.bly <- ifelse( is.na(court.outcome.p.combined.bly), 0,court.outcome.p.combined.bly)
	#for combining regions
	court.outcome.not.p.combined.bly <- ifelse(	court.outcome.p.combined.bly==0,1,0)
	court.outcome.not.gridlock.bly <- ifelse(	court.outcome.gridlock.bly==0,1,0)
	court.outcome.not.flip.bly <- ifelse(court.outcome.flip.bly==0,1,0)
	#pos.taking MODELS
	pos.taking.gridlock.bly <- ifelse(presreg.pos.taking.bly=="gridlock", 1,0)
	pos.taking.flip.bly <- ifelse(presreg.pos.taking.smaller.shift.type.bly== "flip", 1,0)
	pos.taking.flip.bly <- ifelse( is.na(pos.taking.flip.bly ), 0,pos.taking.flip.bly )
	pos.taking.shift.bly <- ifelse(presreg.pos.taking.smaller.shift.type.bly== "shift", 1,0)
	pos.taking.shift.bly <- ifelse( is.na(pos.taking.shift.bly ), 0,pos.taking.shift.bly )
	pos.taking.not.gridlock.bly <- ifelse(	pos.taking.gridlock.bly==0,1,0)
	pos.taking.not.flip.bly <- ifelse(	pos.taking.flip.bly==0,1,0)
	pos.taking.not.shift.bly <- ifelse(	pos.taking.shift.bly==0,1,0)
	#Senate indifference POINT: this applies across models
	sen.ind.point.bly <- 2*sen.median.sim.bly-j5old.sim.bly	 
})

#ncourt.outcome model#
#Model 1: just theory predictions
court.outcome.model.1.bly.intercept <- rep(NA, n.sims)
court.outcome.model.1.bly.gridlock <- rep(NA, n.sims)
court.outcome.model.1.bly.p.predicted <- rep(NA, n.sims)
court.outcome.model.1.bly.flip <- rep(NA, n.sims)
court.outcome.model.1.bly.rsquared <- rep(NA, n.sims)

for (i in 1:n.sims){
	ok <- bailey.nomd$simulation==i
  foo <-	lm(nom.proj.sim.bly ~  	court.outcome.gridlock.bly:j5old.sim.bly + court.outcome.p.combined.bly:pres.sim.bly +	court.outcome.flip.bly:sen.ind.point.bly ,	
		data=bailey.nomd,subset=ok)
	sim.foo <- sim(foo,n.sims=1)
	court.outcome.model.1.bly.intercept[i] <- coef(sim.foo)[1]
	court.outcome.model.1.bly.gridlock[i] <- coef(sim.foo)[2]
	court.outcome.model.1.bly.p.predicted[i] <- coef(sim.foo)[3]
	court.outcome.model.1.bly.flip[i] <- coef(sim.foo)[4]
	court.outcome.model.1.bly.rsquared[i] <- summary(foo)$r.squared
}		

#MODEL 2: now in non-theory predictions
court.outcome.model.2.bly.intercept <- rep(NA, n.sims)
court.outcome.model.2.bly.not.gridlock <- rep(NA, n.sims)
court.outcome.model.2.bly.not.p.predicted <- rep(NA, n.sims)
court.outcome.model.2.bly.not.flip <- rep(NA, n.sims)
court.outcome.model.2.bly.rsquared <- rep(NA, n.sims)

for (i in 1:n.sims){
	ok <- bailey.nomd$simulation==i
  foo <-	lm(nom.proj.sim.bly ~  	
		+ 	court.outcome.not.gridlock.bly:j5old.sim.bly + court.outcome.not.p.combined.bly:pres.sim.bly +	court.outcome.not.flip.bly:sen.ind.point.bly
		, data=bailey.nomd,subset=ok)
	sim.foo <- sim(foo,n.sims=1)
	court.outcome.model.2.bly.intercept[i]    <- coef(sim.foo)[1]
	court.outcome.model.2.bly.not.gridlock[i]    <- coef(sim.foo)[2]
	court.outcome.model.2.bly.not.p.predicted[i] <- coef(sim.foo)[3]
	court.outcome.model.2.bly.not.flip[i]        <- coef(sim.foo)[4]	
	court.outcome.model.2.bly.rsquared[i] <- summary(foo)$r.squared
}		
#########
######################
#pos.taking model
######################
######################
#1) just theory predictions
pos.taking.model.1.bly.intercept <- rep(NA, n.sims)
pos.taking.model.1.bly.gridlock <- rep(NA, n.sims)
pos.taking.model.1.bly.shift <- rep(NA, n.sims)
pos.taking.model.1.bly.flip <- rep(NA, n.sims)
pos.taking.model.1.bly.rsquared <- rep(NA, n.sims)
	
for (i in 1:n.sims){
	ok <- bailey.nomd$simulation==i
	#pos.taking
 foo <-	lm(nom.proj.sim.bly ~  	pos.taking.gridlock.bly:j5old.sim.bly + pos.taking.shift.bly:pres.sim.bly + pos.taking.flip.bly:sen.ind.point.bly,
		 	data=bailey.nomd,subset=ok)
	sim.foo <- sim(foo,n.sims=1)
	pos.taking.model.1.bly.intercept[i] <- coef(sim.foo)[1]
	pos.taking.model.1.bly.gridlock[i] <- coef(sim.foo)[2]
	pos.taking.model.1.bly.shift[i] <- coef(sim.foo)[3]
	pos.taking.model.1.bly.flip[i] <- coef(sim.foo)[4]
	pos.taking.model.1.bly.rsquared[i] <- summary(foo)$r.squared
}		

#2) now non-theory specific predictgions
pos.taking.model.2.bly.intercept <- rep(NA, n.sims)
pos.taking.model.2.bly.not.gridlock <- rep(NA, n.sims)
pos.taking.model.2.bly.not.shift <- rep(NA, n.sims)
pos.taking.model.2.bly.not.flip <- rep(NA, n.sims)
pos.taking.model.2.bly.rsquared <- rep(NA, n.sims)

for (i in 1:n.sims){
	ok <- bailey.nomd$simulation==i
	#pos.taking
 foo <-	lm(nom.proj.sim.bly ~  	pos.taking.not.gridlock.bly:j5old.sim.bly + pos.taking.not.shift.bly:pres.sim.bly + pos.taking.not.flip.bly:sen.ind.point.bly, data=bailey.nomd,subset=ok)
	sim.foo <- sim(foo,n.sims=1)
	pos.taking.model.2.bly.intercept[i] <- coef(foo)[1]
	pos.taking.model.2.bly.not.gridlock[i] <-  coef(foo)[2]
	pos.taking.model.2.bly.not.shift[i] <- coef(foo)[3]
  pos.taking.model.2.bly.not.flip[i]  <-	 coef(foo)[4]
	pos.taking.model.2.bly.rsquared[i] <- summary(foo)$r.squared
}
#REPORT RESULTS FOR TABLE 3
#NEARLY court.outcome
	#CORRECT MODELS
#MODEL 1 (DW)
quantile(court.outcome.model.1.dw.intercept , c(lo,.5,hi))
quantile(court.outcome.model.1.dw.gridlock, c(lo,.5,hi))
quantile(	court.outcome.model.1.dw.p.predicted, c(lo,.5,hi))
quantile(	court.outcome.model.1.dw.flip , c(lo,.5,hi), na.rm=T)
mean(	court.outcome.model.1.dw.rsquared)
#MODEL 2 (BLY)
quantile(court.outcome.model.1.bly.intercept , c(lo,.5,hi))
quantile(court.outcome.model.1.bly.gridlock, c(lo,.5,hi))
quantile(	court.outcome.model.1.bly.p.predicted, c(lo,.5,hi))
quantile(	court.outcome.model.1.bly.flip , c(lo,.5,hi), na.rm=T)
mean(	court.outcome.model.1.bly.rsquared)
#PLACEBO MODELS
#Model 3 (DW)
quantile(court.outcome.model.2.dw.intercept , c(lo,.5,hi))
quantile(court.outcome.model.2.dw.not.gridlock, c(lo,.5,hi))
quantile(court.outcome.model.2.dw.not.p.predicted, c(lo,.5,hi))
quantile(court.outcome.model.2.dw.not.flip , c(lo,.5,hi), na.rm=T)
mean(	court.outcome.model.2.dw.rsquared)
#Model 4 (BLY)
quantile(court.outcome.model.2.bly.intercept , c(lo,.5,hi))
quantile(court.outcome.model.2.bly.not.gridlock, c(lo,.5,hi))
quantile(court.outcome.model.2.bly.not.p.predicted, c(lo,.5,hi))
quantile(court.outcome.model.2.bly.not.flip , c(lo,.5,hi), na.rm=T)
mean(	court.outcome.model.2.bly.rsquared)
mean(	court.outcome.model.1.bly.p.predicted>court.outcome.model.2.bly.not.p.predicted)

#pos.taking
#CORRECT
#Model 5 (DW)
quantile(	pos.taking.model.1.dw.intercept, c(lo,.5,hi))
quantile(	pos.taking.model.1.dw.gridlock, c(lo,.5,hi))
quantile(	pos.taking.model.1.dw.shift, c(lo,.5,hi))
quantile(	pos.taking.model.1.dw.flip, c(lo,.5,hi))
mean(	pos.taking.model.1.dw.rsquared)
#Model 6 (Bailey)
quantile(	pos.taking.model.1.bly.intercept, c(lo,.5,hi))
quantile(	pos.taking.model.1.bly.gridlock, c(lo,.5,hi))
quantile(	pos.taking.model.1.bly.shift, c(lo,.5,hi))
quantile(	pos.taking.model.1.bly.flip, c(lo,.5,hi))
mean(	pos.taking.model.1.bly.rsquared)
#Model (7) 
quantile( pos.taking.model.2.dw.intercept, c(lo,.5,hi))
quantile(	pos.taking.model.2.dw.not.gridlock, c(lo,.5,hi))
quantile(	pos.taking.model.2.dw.not.shift, c(lo,.5,hi))
quantile(	pos.taking.model.2.dw.not.flip, c(lo,.5,hi))
mean(	pos.taking.model.2.dw.rsquared)
#Model (8)
quantile( pos.taking.model.2.bly.intercept, c(lo,.5,hi))
quantile(	pos.taking.model.2.bly.not.gridlock, c(lo,.5,hi))
quantile(	pos.taking.model.2.bly.not.shift, c(lo,.5,hi))
quantile(	pos.taking.model.2.bly.not.flip, c(lo,.5,hi))
mean(	pos.taking.model.2.bly.rsquared)

######################
#ANALYSES FOR APPENDIX
####################
#Section A.2--mapping justices into DW-nominate
#take means across sims to ignore uncertainty
pres.ideal <- with(nomd, tapply(arctan.pres.sim.dw, nominee,mean))
first.MQ <- with(nomd, tapply(simulated.1st.MQ, nominee,mean))
nom.temp <- names(pres.ideal)
foo <- data.frame(pres.ideal, first.MQ, nom.temp)
#now run regressoins
#1st, unconstrained nominees from Moraski and Shipan
keep <- c("warren", "brennan", "whittaker", "stewart", "white", "goldberg", "marshall", "burger", "haysnworth", "carswell", "blackmun",
	"powell", "rehnquist1", "rehnquist2", "scalia", "ginsburg_r", "breyer")
model.A1.1  <- lm(pres.ideal ~first.MQ, data=foo, subset=nom.temp %in% keep)
arm::display(model.A1.1)
#2nd, all nominees from Moraski
keep <- c("warren", "brennan", "whittaker", "stewart", "white", "goldberg", "marshall", "burger", "haysnworth", "carswell", "blackmun",
	"powell", "rehnquist1", "rehnquist2", "scalia", "ginsburg_r", "breyer", "stevens","bork", "kennedy", "clark", "minton", "harlan", 
	"fortas1", "fortas2", "souter", "thomas")
model.A1.2  <- lm(pres.ideal ~first.MQ, data=foo, subset=nom.temp %in% keep)
arm::display(model.A1.2)
#3rd, all nominees in our sample
model.A1.3  <- lm(pres.ideal ~first.MQ, data=foo)
arm::display(model.A1.3)
apsrtable::apsrtable(model.A1.1,model.A1.2, model.A1.3)



#REPLICATE RELEVANT ANALYSES USING FILIBUSTER PIVOT (SECTION A.5)
###################################################
#2) "Off-sides/AGGRESIVE MISTAKE" moaving new median/nominee. farther away from sen. than old median
#court.outcome AND pos.taking--DEFINE IN TERMS OF BOTH NEW MEDIAN AND OLD MEDIAN
###################################################
#DW
#1st, in terms of new median (court.outcome)
nomd.pivot <- nomd #create new DATAFRAMES to distinguish analyses (mainly regressions)
bailey.nomd.pivot <- bailey.nomd
nom.level.stats.pivot <- nom.level.stats
bailey.nom.level.stats.pivot <- bailey.nom.level.stats

nomd.pivot$agg.mistake.court.outcome.amount.dw <- abs( nomd.pivot$sen.pivot.sim.dw- nomd.pivot$j5new.sim.dw) - abs(  nomd.pivot$sen.pivot.sim.dw- nomd.pivot$j5old.sim.dw) 
nom.level.stats.pivot$agg.mistake.court.outcome.distance.lo.dw <- tapply(nomd.pivot$agg.mistake.court.outcome.amount.dw, nomd.pivot$nominee, quantile,probs=c(lo.ci))
nom.level.stats.pivot$agg.mistake.court.outcome.distance.med.dw <-  tapply(nomd.pivot$agg.mistake.court.outcome.amount.dw, nomd.pivot$nominee, quantile,probs=c(.5))
nom.level.stats.pivot$agg.mistake.court.outcome.distance.hi.dw <- tapply(nomd.pivot$agg.mistake.court.outcome.amount.dw, nomd.pivot$nominee, quantile,probs=c(hi.ci))
nom.level.stats.pivot$agg.mistake.court.outcome.mean.prob.dw <- tapply(nomd.pivot$agg.mistake.court.outcome.amount.dw > 0, nomd.pivot$nominee,mean)#prob of an own goal

y.axis.text <- .7

#which nominees were confirmed among those with a 95\% probability of being an aggressive mistake
with(nom.level.stats.pivot, table(confirmed[agg.mistake.court.outcome.mean.prob.dw>=.95]))
#DW
pdf("Figure_A5_bottom_left.pdf", height=7, width=5.5)

y <- 1:length(nom.level.stats.pivot$agg.mistake.court.outcome.mean.prob.dw )
x <- sort (nom.level.stats.pivot$agg.mistake.court.outcome.distance.med.dw)
x.lo <- nom.level.stats.pivot$agg.mistake.court.outcome.distance.lo.dw[order(nom.level.stats.pivot$agg.mistake.court.outcome.distance.med.dw)]
x.hi <- nom.level.stats.pivot$agg.mistake.court.outcome.distance.hi.dw[order(nom.level.stats.pivot$agg.mistake.court.outcome.distance.med.dw)]
prob.agg.mistake.court.outcome.sorted <- nom.level.stats.pivot$agg.mistake.court.outcome.mean.prob.dw [order(nom.level.stats.pivot$agg.mistake.court.outcome.distance.med.dw, decreasing=F)]
nominee.sorted <- nominee.title.dw[order(nom.level.stats.pivot$agg.mistake.court.outcome.distance.med.dw)]
confirmed.sorted <- nom.level.stats.pivot$confirmed[order(nom.level.stats.pivot$agg.mistake.court.outcome.distance.med.dw)]
prob.text <- .8

par(mfrow=c(1,1),mar=c(4,4.5,0,0))
plot(x,y , xlim=c(-.8,.8), pch=19, col="black", xlab="", ylab="",axes=F, type="n")

polygon(x=c(0,0,par()$usr[2], par()$usr[2]),
             y=par()$usr[c(3,4,4,3)],  ## vertical limits of plot
             col=gray(shading),   ## desired color
             border=F)         ## no border


segments(x.lo, y, x.hi, y)
points(x[confirmed.sorted==1], y[confirmed.sorted==1], pch=19)
points(x[confirmed.sorted==0], y[confirmed.sorted==0], pch=21, bg="white")

axis(1, at = seq(-.8,.8,.2), 
     labels=c("-.8","-.6",  "-.4", "-.2", 0, ".2", ".4", ".6",".8"), mgp=c(2,.5,0), cex.axis=x.axis)
axis(2, at =y , label=nominee.sorted, mgp=c(2,.5,0), cex.axis=y.axis.text, las=2)
abline(v=0, col="gray90")
#add probability for nominees greater than zero
ok <- x>0
text(x.hi[ok]+.08,y[ok], round(prob.agg.mistake.court.outcome.sorted[ok],2), col="red", cex=prob.text)
text(.5, 5.1, "Aggressive\nmistakes,\nNOMINATE\n(in terms of\nnew median)", cex=1.3)

#mtext("| Distance from new median to senate median | -\n| distance from old median to senate median | ", 1, line=2.5)#, adj=1)
text(0,-5,expression(paste("| ", italic(s[fp])- italic(j)[5]^1 , " | - | ", italic(s[fp])- italic(j)[5]^0, " |"   )), cex=1.3, xpd=NA)

dev.off()

#2nd, agg.mistake in terms of nominee (pos.taking)
#DW
nomd.pivot$agg.mistake.pos.taking.amount.dw <- abs( nomd.pivot$sen.pivot.sim.dw- nomd.pivot$nom.proj.sim.dw) - abs(nomd.pivot$sen.pivot.sim.dw- nomd.pivot$j5old.sim.dw) 
nom.level.stats.pivot$agg.mistake.pos.taking.distance.lo.dw <- tapply(nomd.pivot$agg.mistake.pos.taking.amount.dw, nomd.pivot$nominee, quantile,probs=c(lo.ci))
nom.level.stats.pivot$agg.mistake.pos.taking.distance.med.dw <-  tapply(nomd.pivot$agg.mistake.pos.taking.amount.dw, nomd.pivot$nominee, quantile,probs=c(.5))
nom.level.stats.pivot$agg.mistake.pos.taking.distance.hi.dw <- tapply(nomd.pivot$agg.mistake.pos.taking.amount.dw, nomd.pivot$nominee, quantile,probs=c(hi.ci))
nom.level.stats.pivot$agg.mistake.pos.taking.mean.prob.dw <- tapply(nomd.pivot$agg.mistake.pos.taking.amount.dw > 0, nomd.pivot$nominee,mean)#prob of an own goal

pdf("Figure_A5_top_left.pdf", height=7, width=5.5)

y <- 1:length(nom.level.stats.pivot$agg.mistake.pos.taking.mean.prob.dw )
x <- sort (nom.level.stats.pivot$agg.mistake.pos.taking.distance.med.dw)
x.lo <- nom.level.stats.pivot$agg.mistake.pos.taking.distance.lo.dw[order(nom.level.stats.pivot$agg.mistake.pos.taking.distance.med.dw)]
x.hi <- nom.level.stats.pivot$agg.mistake.pos.taking.distance.hi.dw[order(nom.level.stats.pivot$agg.mistake.pos.taking.distance.med.dw)]
prob.agg.mistake.pos.taking.sorted <- nom.level.stats.pivot$agg.mistake.pos.taking.mean.prob.dw [order(nom.level.stats.pivot$agg.mistake.pos.taking.distance.med.dw, decreasing=F)]
nominee.sorted <- nominee.title.dw[order(nom.level.stats.pivot$agg.mistake.pos.taking.distance.med.dw)]
confirmed.sorted <- nom.level.stats.pivot$confirmed[order(nom.level.stats.pivot$agg.mistake.pos.taking.distance.med.dw)]

par(mfrow=c(1,1),mar=c(4,4.5,0,0))
plot(x,y , xlim=c(-.8,.8), pch=19, col="black", xlab="", ylab="",axes=F, type="n")
polygon(x=c(0,0,par()$usr[2], par()$usr[2]),
             y=par()$usr[c(3,4,4,3)],  ## vertical limits of plot
             col=gray(shading),   ## desired color
             border=F)         ## no border

segments(x.lo, y, x.hi, y)
points(x[confirmed.sorted==1], y[confirmed.sorted==1], pch=19)
points(x[confirmed.sorted==0], y[confirmed.sorted==0], pch=21, bg="white")

axis(1, at = seq(-.8,.8,.2), 
     labels=c("-.8","-.6",  "-.4", "-.2", 0, ".2", ".4", ".6",".8"), mgp=c(2,.5,0), cex.axis=x.axis)
axis(2, at =y , label=nominee.sorted, mgp=c(2,.5,0), cex.axis=y.axis.text, las=2)
abline(v=0, col="gray90")
#mtext("| Distance from Nominee to senate median | -\n| distance from old median to senate median | ", 1, line=2.5)#, adj=1)
#add probability for nominees greater than zero
ok <- x>0
text(x.hi[ok]+.1,y[ok], round(prob.agg.mistake.pos.taking.sorted[ok],2), col="red")
text(.5, 4.4, "Aggressive\nmistakes,\nNOMINATE\n(in terms\nof nominee)", cex=1.3)
text(0,-5,expression(paste("| ",  italic(s[fp]) - italic(n), " | - | ", italic(s[fp])- italic(j)[5]^0 , " |"   )), cex=1.3, xpd=NA)

dev.off()

#bly
#1st, in terms of new median (court.outcome)
bailey.nomd.pivot$agg.mistake.court.outcome.amount.bly <- abs( bailey.nomd.pivot$sen.pivot.sim.bly- bailey.nomd.pivot$j5new.sim.bly) - abs(  bailey.nomd.pivot$sen.pivot.sim.bly- bailey.nomd.pivot$j5old.sim.bly) 
bailey.nom.level.stats.pivot$agg.mistake.court.outcome.distance.lo.bly <- tapply(bailey.nomd.pivot$agg.mistake.court.outcome.amount.bly, bailey.nomd.pivot$nominee, quantile,probs=c(lo.ci))
bailey.nom.level.stats.pivot$agg.mistake.court.outcome.distance.med.bly <-  tapply(bailey.nomd.pivot$agg.mistake.court.outcome.amount.bly, bailey.nomd.pivot$nominee, quantile,probs=c(.5))
bailey.nom.level.stats.pivot$agg.mistake.court.outcome.distance.hi.bly <- tapply(bailey.nomd.pivot$agg.mistake.court.outcome.amount.bly, bailey.nomd.pivot$nominee, quantile,probs=c(hi.ci))
bailey.nom.level.stats.pivot$agg.mistake.court.outcome.mean.prob.bly <- tapply(bailey.nomd.pivot$agg.mistake.court.outcome.amount.bly > 0, bailey.nomd.pivot$nominee,mean)#prob of an own goal

#which nominees were confirmed among those with a 95\% probability of being an aggressive mistake
with(bailey.nom.level.stats.pivot, table(confirmed[agg.mistake.court.outcome.mean.prob.bly>=.95]))

pdf("Figure_A5_bottom_right.pdf", height=7, width=5.5)

y <- 1:length(bailey.nom.level.stats.pivot$agg.mistake.court.outcome.mean.prob.bly )
x <- sort (bailey.nom.level.stats.pivot$agg.mistake.court.outcome.distance.med.bly)
x.lo <- bailey.nom.level.stats.pivot$agg.mistake.court.outcome.distance.lo.bly[order(bailey.nom.level.stats.pivot$agg.mistake.court.outcome.distance.med.bly)]
x.hi <- bailey.nom.level.stats.pivot$agg.mistake.court.outcome.distance.hi.bly[order(bailey.nom.level.stats.pivot$agg.mistake.court.outcome.distance.med.bly)]
prob.agg.mistake.court.outcome.sorted <- bailey.nom.level.stats.pivot$agg.mistake.court.outcome.mean.prob.bly [order(bailey.nom.level.stats.pivot$agg.mistake.court.outcome.distance.med.bly, decreasing=F)]
nominee.sorted <- nominee.title.bly[order(bailey.nom.level.stats.pivot$agg.mistake.court.outcome.distance.med.bly)]
confirmed.sorted <- bailey.nom.level.stats.pivot$confirmed[order(bailey.nom.level.stats.pivot$agg.mistake.court.outcome.distance.med.bly)]
prob.text <- .8

par(mfrow=c(1,1),mar=c(4,4.5,0,0))
plot(x,y , xlim=c(-1,1.4), pch=19, col="black", xlab="", ylab="",axes=F, type="n")
polygon(x=c(0,0,par()$usr[2], par()$usr[2]),
             y=par()$usr[c(3,4,4,3)],  ## vertical limits of plot
             col=gray(shading),   ## desired color
             border=F)         ## no border

segments(x.lo, y, x.hi, y)
points(x[confirmed.sorted==1], y[confirmed.sorted==1], pch=19)
points(x[confirmed.sorted==0], y[confirmed.sorted==0], pch=21, bg="white")

axis(1, at = seq(-1,1.4,.2), 
     labels=c("-1","-.8","-.6",  "-.4", "-.2", 0, ".2", ".4", ".6",".8","1","1.2","1.4"), mgp=c(2,.5,0), cex.axis=x.axis)
axis(2, at =y , label=nominee.sorted, mgp=c(2,.5,0), cex.axis=y.axis.text, las=2)
abline(v=0, col="gray90")
#mtext("| Distance from new median to senate median | -\n| distance from old median to senate median | ", 1, line=2.5)#, adj=1)
#add probability for nominees greater than zero
ok <- x>0
text(x.hi[ok]+.08,y[ok], round(prob.agg.mistake.court.outcome.sorted[ok],2), col="red", cex=prob.text)
text(.8, 3.5, "Aggressive\nmistakes, Bailey \n(in terms of\nnew median)", cex=1.3)
text(0,-3,expression(paste("| ", italic(s[fp])- italic(j)[5]^1 , " | - | ", italic(s[fp])- italic(j)[5]^0, " |"   )), cex=1.3, xpd=NA)

dev.off()

#now in terms of nominee
bailey.nomd.pivot$agg.mistake.pos.taking.amount.bly <- abs( bailey.nomd.pivot$sen.pivot.sim.bly- bailey.nomd.pivot$nom.proj.sim.bly) - abs(bailey.nomd.pivot$sen.pivot.sim.bly- bailey.nomd.pivot$j5old.sim.bly) 
bailey.nom.level.stats.pivot$agg.mistake.pos.taking.distance.lo.bly <- tapply(bailey.nomd.pivot$agg.mistake.pos.taking.amount.bly, bailey.nomd.pivot$nominee, quantile,probs=c(lo.ci))
bailey.nom.level.stats.pivot$agg.mistake.pos.taking.distance.med.bly <-  tapply(bailey.nomd.pivot$agg.mistake.pos.taking.amount.bly, bailey.nomd.pivot$nominee, quantile,probs=c(.5))
bailey.nom.level.stats.pivot$agg.mistake.pos.taking.distance.hi.bly <- tapply(bailey.nomd.pivot$agg.mistake.pos.taking.amount.bly, bailey.nomd.pivot$nominee, quantile,probs=c(hi.ci))
bailey.nom.level.stats.pivot$agg.mistake.pos.taking.mean.prob.bly <- tapply(bailey.nomd.pivot$agg.mistake.pos.taking.amount.bly > 0, bailey.nomd.pivot$nominee,mean)#prob of an own goal

pdf("Figure_A5_top_right.pdf", height=7, width=5.5)

y <- 1:length(bailey.nom.level.stats.pivot$agg.mistake.pos.taking.mean.prob.bly )
x <- sort (bailey.nom.level.stats.pivot$agg.mistake.pos.taking.distance.med.bly)
x.lo <- bailey.nom.level.stats.pivot$agg.mistake.pos.taking.distance.lo.bly[order(bailey.nom.level.stats.pivot$agg.mistake.pos.taking.distance.med.bly)]
x.hi <- bailey.nom.level.stats.pivot$agg.mistake.pos.taking.distance.hi.bly[order(bailey.nom.level.stats.pivot$agg.mistake.pos.taking.distance.med.bly)]
prob.agg.mistake.pos.taking.sorted <- bailey.nom.level.stats.pivot$agg.mistake.pos.taking.mean.prob.bly [order(bailey.nom.level.stats.pivot$agg.mistake.pos.taking.distance.med.bly, decreasing=F)]
nominee.sorted <- nominee.title.bly[order(bailey.nom.level.stats.pivot$agg.mistake.pos.taking.distance.med.bly)]
confirmed.sorted <- bailey.nom.level.stats.pivot$confirmed[order(bailey.nom.level.stats.pivot$agg.mistake.pos.taking.distance.med.bly)]

par(mfrow=c(1,1),mar=c(4,4.5,0,0))
plot(x,y , xlim=c(-1,1.4), pch=19, col="black", xlab="", ylab="",axes=F, type="n")
polygon(x=c(0,0,par()$usr[2], par()$usr[2]),
             y=par()$usr[c(3,4,4,3)],  ## vertical limits of plot
             col=gray(shading),   ## desired color
             border=F)         ## no border
segments(x.lo, y, x.hi, y)
points(x[confirmed.sorted==1], y[confirmed.sorted==1], pch=19)
points(x[confirmed.sorted==0], y[confirmed.sorted==0], pch=21, bg="white")

axis(1, at = seq(-1,1.4,.2), 
     labels=c("-1","-.8","-.6",  "-.4", "-.2", 0, ".2", ".4", ".6",".8","1","1.2","1.4"), mgp=c(2,.5,0), cex.axis=x.axis)
axis(2, at =y , label=nominee.sorted, mgp=c(2,.5,0), cex.axis=y.axis.text, las=2)
abline(v=0, col="gray90")
#mtext("| Distance from Nominee to senate median | -\n| distance from old median to senate median | ", 1, line=2.5)#, adj=1)
#add probability for nominees greater than zero
ok <- x>0
text(x.hi[ok]+.1,y[ok], round(prob.agg.mistake.pos.taking.sorted[ok],2), col="red")
text(.8, 3.5, "Aggressive\nmistakes, Bailey\n(in terms\nof nominee)", cex=1.3)
text(0,-3,expression(paste("| ",  italic(s[fp]) - italic(n), " | - | ", italic(s[fp])- italic(j)[5]^0 , " |"   )), cex=1.3, xpd=NA)

dev.off()

##########################################################################################
#now calculate predicted location of nominee, calculating regime in each sim#
#first, get SENATE REGIMES
#make senator cutpoints to divide A/B and C/D
#DW
nomd.pivot$j4j5mid.sim.dw <- (nomd.pivot$j4old.sim.dw + nomd.pivot$j5old.sim.dw) / 2
nomd.pivot$j5j6mid.sim.dw <- (nomd.pivot$j5old.sim.dw + nomd.pivot$j6old.sim.dw) / 2 

nomd.pivot <- within(nomd.pivot, {
	sen.regime.dw.pivot <- ifelse( sen.pivot.sim.dw  <  j4j5mid.sim.dw, "A", 
								ifelse( (j4j5mid.sim.dw <= sen.pivot.sim.dw) & (sen.pivot.sim.dw < j5old.sim.dw), "B",
								ifelse( (j5old.sim.dw <=  sen.pivot.sim.dw)  & (sen.pivot.sim.dw <=  j5j6mid.sim.dw), "C",
								ifelse( sen.pivot.sim.dw > j5j6mid.sim.dw, "D", NA))))
	})	

####NOW SET UP DIRECT TESTS (Including 3rd and 4th robust regimes)
#FIRST, FOR court.outcome AND NEARLY court.outcome, DESCRIBE REGIMES
nomd.pivot <- within(nomd.pivot, {
	presreg.court.outcome.dw.pivot <- ifelse( above.dw==1 &  j.exit.sim.dw > j5old.sim.dw,   "restoring",
									ifelse( ( above.dw==1 &  j.exit.sim.dw <= j5old.sim.dw & (sen.regime.dw.pivot == "A" | sen.regime.dw.pivot=="B")),  "gridlock",
									ifelse( above.dw== 1 & ( j.exit.sim.dw  <= j5old.sim.dw) & ( sen.regime.dw.pivot  == "C"), "smaller.shift",
									ifelse(above.dw==1 & ( j.exit.sim.dw  <= j5old.sim.dw) & ( sen.regime.dw.pivot == "D"), "maxshift",
									ifelse( above.dw==0 &  j.exit.sim.dw < j5old.sim.dw,   "restoring",
									ifelse( ( above.dw==0 &  j.exit.sim.dw >= j5old.sim.dw & (sen.regime.dw.pivot == "C" | sen.regime.dw.pivot=="D")),  "gridlock",
									ifelse( above.dw== 0 & ( j.exit.sim.dw  >= j5old.sim.dw) & ( sen.regime.dw.pivot  == "B"), "smaller.shift",
									ifelse(above.dw==0 & ( j.exit.sim.dw >= j5old.sim.dw) & ( sen.regime.dw.pivot == "A"), "maxshift",
									NA))))))))
#now get point prediction for nearly court.outcome
ncourt.outcome.point.pred.dw.pivot <- ifelse(	presreg.court.outcome.dw.pivot=="restoring", pres.sim.dw,
										ifelse(	presreg.court.outcome.dw.pivot=="gridlock", j5old.sim.dw,
										ifelse(	presreg.court.outcome.dw.pivot=="smaller.shift" & above.dw==1, apply(  matrix(c(pres.sim.dw, (2*sen.pivot.sim.dw-j5old.sim.dw)),ncol=2) , 1, min),						
										ifelse(	presreg.court.outcome.dw.pivot=="smaller.shift" & above.dw==0, apply(  matrix(c(pres.sim.dw, (2*sen.pivot.sim.dw-j5old.sim.dw)),ncol=2) , 1, max),
										ifelse(	presreg.court.outcome.dw.pivot=="maxshift",pres.sim.dw, NA)))))			
									
# now define whether nominee is flip or shift within "smaller.shift"
# shift is if president is interior" to 2S-j -- can get his ideal point  in nearly court.outcome and mixed game
# flip  is if president is exterior to to 2S-j -- can only get S's reflection point  in nearly court.outcome and mixed game 3
#there are only about 300 simulations where shift is met						
	presreg.court.outcome.smaller.shift.type.dw.pivot <- 
			ifelse(above.dw==1 & 	presreg.court.outcome.dw.pivot=="smaller.shift"& pres.sim.dw <= 2*sen.pivot.sim.dw-j5old.sim.dw, "shift",
			ifelse(above.dw==1 & 	presreg.court.outcome.dw.pivot=="smaller.shift"& pres.sim.dw > 2*sen.pivot.sim.dw-j5old.sim.dw, "flip",
			ifelse(above.dw==0 & 	presreg.court.outcome.dw.pivot=="smaller.shift"& pres.sim.dw >= 2*sen.pivot.sim.dw-j5old.sim.dw, "shift",
			ifelse(above.dw==0 & 	presreg.court.outcome.dw.pivot=="smaller.shift"& pres.sim.dw < 2*sen.pivot.sim.dw-j5old.sim.dw, "flip", NA))))	
#NOW DO pos.taking REGIMES
		presreg.pos.taking.dw.pivot <- ifelse( above.dw==1 & (sen.regime.dw.pivot == "A" | sen.regime.dw.pivot=="B"), "gridlock",
									  	ifelse( above.dw==1 & (sen.regime.dw.pivot == "C" | sen.regime.dw.pivot=="D"), "smaller.shift",
										  ifelse( above.dw==0 & (sen.regime.dw.pivot == "C" | sen.regime.dw.pivot=="D"), "gridlock",
										  ifelse( above.dw==0 & (sen.regime.dw.pivot == "A" | sen.regime.dw.pivot=="B"), "smaller.shift", NA ))))
#get point prediction for pos.taking
   pos.taking.point.pred.dw.pivot <- ifelse(presreg.pos.taking.dw.pivot=="gridlock", j5old.sim.dw,
										ifelse(	presreg.pos.taking.dw.pivot=="smaller.shift" & above.dw==1, apply(  matrix(c(pres.sim.dw, (2*sen.pivot.sim.dw-j5old.sim.dw)),ncol=2) , 1, min),						
										ifelse(	presreg.pos.taking.dw.pivot=="smaller.shift" & above.dw==0, apply(  matrix(c(pres.sim.dw, (2*sen.pivot.sim.dw-j5old.sim.dw)),ncol=2) , 1, max), NA)))																 
 	 presreg.pos.taking.smaller.shift.type.dw.pivot <- 
			ifelse(above.dw==1 & 	presreg.pos.taking.dw.pivot=="smaller.shift" & pres.sim.dw <= 2*sen.pivot.sim.dw-j5old.sim.dw, "shift",
			ifelse(above.dw==1 & 	presreg.pos.taking.dw.pivot=="smaller.shift" & pres.sim.dw > 2*sen.pivot.sim.dw-j5old.sim.dw, "flip",
			ifelse(above.dw==0 & 	presreg.pos.taking.dw.pivot=="smaller.shift" & pres.sim.dw >= 2*sen.pivot.sim.dw-j5old.sim.dw, "shift",
			ifelse(above.dw==0 & 	presreg.pos.taking.dw.pivot=="smaller.shift" & pres.sim.dw < 2*sen.pivot.sim.dw-j5old.sim.dw, "flip", NA))))	
	})

table(nomd.pivot$sen.regime.dw.pivot)/length(nomd.pivot$sen.regime.dw)
table(nomd$sen.regime.dw)/length(nomd.pivot$sen.regime.dw)
mean(nomd.pivot$sen.regime.dw)

tapply(nomd.pivot$sen.regime.dw, nomd.pivot$nominee, table)
tapply(c(nomd.pivot$presreg.court.outcome.dw,nomd.pivot$sen.regime.dw), c(nomd.pivot$nominee,nomd.pivot$nominee) , table)
tapply(c(nomd.pivot$presreg.pos.taking.dw,nomd.pivot$sen.regime.dw), c(nomd.pivot$nominee,nomd.pivot$nominee) , table)

#make matrix with frequencies of all categories by nominee
library(zoo)
court.outcome.regime.dist.dw.pivot  <- data.frame(t(do.call(merge, lapply(tapply(nomd.pivot$presreg.court.outcome.dw.pivot, nomd.pivot$nominee, table), function(x) zoo(c(x), rownames(x))))))
court.outcome.regime.dist.dw.pivot <- within(court.outcome.regime.dist.dw.pivot, {
	nominee <- rownames(court.outcome.regime.dist.dw.pivot)
	court.outcome.regime.modal.dw.pivot  <- apply(court.outcome.regime.dist.dw.pivot ,1,which.max)
	court.outcome.regime.modal.dw.pivot <- ifelse(court.outcome.regime.modal.dw.pivot==1, "flip", ifelse(court.outcome.regime.modal.dw.pivot==2, "gridlock", 
											 ifelse(court.outcome.regime.modal.dw.pivot=="3", "maxshift", ifelse(court.outcome.regime.modal.dw.pivot==4,"restoring", NA))))
  court.outcome.regime.max.percent.dw.pivot <- apply(cbind(smaller.shift, gridlock, maxshift, restoring),1,max, na.rm=T)  / n.sims
})
head(court.outcome.regime.dist.dw.pivot)
hist(court.outcome.regime.dist.dw.pivot$court.outcome.regime.max.percent.dw.pivot)
table(court.outcome.regime.dist.dw.pivot$court.outcome.regime.max.percent.dw.pivot >.9)

#can we test flip or shift prediction?
table(court.outcome.regime.dist.dw.pivot$smaller.shift/n.sims>.5)#NO-

pos.taking.regime.dist.dw.pivot  <- data.frame(t(do.call(merge, lapply(tapply(nomd.pivot$presreg.pos.taking.dw.pivot, nomd.pivot$nominee, table), function(x) zoo(c(x), rownames(x))))))
pos.taking.regime.dist.dw.pivot <- within(pos.taking.regime.dist.dw.pivot, {
  nominee <- rownames(pos.taking.regime.dist.dw.pivot )
	pos.taking.regime.modal.dw.pivot  <- apply(pos.taking.regime.dist.dw.pivot ,1,which.max)
	pos.taking.regime.modal.dw.pivot <- ifelse(pos.taking.regime.modal.dw.pivot==1, "smaller.shift", ifelse(pos.taking.regime.modal.dw.pivot ==2, "gridlock",NA))
  pos.taking.regime.max.percent.dw.pivot <- apply(cbind(smaller.shift, gridlock),1,max, na.rm=T)  / n.sims
})
table(pos.taking.regime.dist.dw.pivot$pos.taking.regime.max.percent.dw.pivot>.8)

#separately, define "lower left" in terms of pres theory figure paper for testing 3rd robust prediction
#need to account for president
nomd.pivot <- within(nomd.pivot, {
	lower.left.dw.pivot <- ifelse( (above.dw ==1  & (sen.regime.dw.pivot == "A" | sen.regime.dw.pivot == "B") &  j.exit.sim.dw <= j5old.sim.dw )  | 
		(above.dw ==0  & (sen.regime.dw.pivot == "C" | sen.regime.dw.pivot == "D") &  j.exit.sim.dw > j5old.sim.dw) ,1,0)
})
prob.lower.left.dw.pivot <-   tapply (nomd.pivot$lower.left.dw, nomd.pivot$nominee,mean, na.omit=F)				

temp.df.pivot <- data.frame(nominee=court.outcome.regime.dist.dw.pivot$nominee,court.outcome.regime.modal.dw.pivot= court.outcome.regime.dist.dw.pivot$court.outcome.regime.modal.dw.pivot,
		court.outcome.regime.max.percent.dw.pivot = court.outcome.regime.dist.dw.pivot$court.outcome.regime.max.percent.dw.pivot ,
	 pos.taking.regime.modal.dw.pivot= pos.taking.regime.dist.dw.pivot$pos.taking.regime.modal.dw.pivot, 
		pos.taking.regime.max.percent.dw.pivot= pos.taking.regime.dist.dw.pivot$pos.taking.regime.max.percent.dw.pivot ,
		prob.lower.left.dw.pivot=prob.lower.left.dw.pivot )
	
#join to nom-data stats
nomd.pivot <- join(nomd.pivot,temp.df.pivot )
table(nomd.pivot$court.outcome.regime.modal.dw.pivot,nomd.pivot$pos.taking.regime.modal.dw.pivot )
head(nomd.pivot)


############
#BAiley
bailey.nomd.pivot$j4j5mid.sim.bly <- (bailey.nomd.pivot$j4old.sim.bly + bailey.nomd.pivot$j5old.sim.bly) / 2
bailey.nomd.pivot$j5j6mid.sim.bly <- (bailey.nomd.pivot$j5old.sim.bly + bailey.nomd.pivot$j6old.sim.bly) / 2 

bailey.nomd.pivot <- within(bailey.nomd.pivot, {
	sen.regime.bly.pivot <- ifelse( sen.pivot.sim.bly  <  j4j5mid.sim.bly, "A", 
								ifelse( (j4j5mid.sim.bly <= sen.pivot.sim.bly) & (sen.pivot.sim.bly < j5old.sim.bly), "B",
								ifelse( (j5old.sim.bly <=  sen.pivot.sim.bly)  & (sen.pivot.sim.bly <=  j5j6mid.sim.bly), "C",
								ifelse( sen.pivot.sim.bly > j5j6mid.sim.bly, "D", NA))))
	})	

####NOW SET UP DIRECT TESTS (Including 3rd and 4th robust regimes)
#FIRST, FOR court.outcome AND NEARLY court.outcome, DESCRIBE REGIMES
bailey.nomd.pivot <- within(bailey.nomd.pivot, {
	presreg.court.outcome.bly.pivot <- ifelse( above.bly==1 &  j.exit.sim.bly > j5old.sim.bly,   "restoring",
									ifelse( ( above.bly==1 &  j.exit.sim.bly <= j5old.sim.bly & (sen.regime.bly.pivot == "A" | sen.regime.bly.pivot=="B")),  "gridlock",
									ifelse( above.bly== 1 & ( j.exit.sim.bly  <= j5old.sim.bly) & ( sen.regime.bly.pivot  == "C"), "smaller.shift",
									ifelse(above.bly==1 & ( j.exit.sim.bly  <= j5old.sim.bly) & ( sen.regime.bly.pivot == "D"), "maxshift",
									ifelse( above.bly==0 &  j.exit.sim.bly < j5old.sim.bly,   "restoring",
									ifelse( ( above.bly==0 &  j.exit.sim.bly >= j5old.sim.bly & (sen.regime.bly.pivot == "C" | sen.regime.bly.pivot=="D")),  "gridlock",
									ifelse( above.bly== 0 & ( j.exit.sim.bly  >= j5old.sim.bly) & ( sen.regime.bly.pivot  == "B"), "smaller.shift",
									ifelse(above.bly==0 & ( j.exit.sim.bly >= j5old.sim.bly) & ( sen.regime.bly.pivot == "A"), "maxshift",
									NA))))))))
#now get point prediction for nearly court.outcome
ncourt.outcome.point.pred.bly.pivot <- ifelse(	presreg.court.outcome.bly.pivot=="restoring", pres.sim.bly,
										ifelse(	presreg.court.outcome.bly.pivot=="gridlock", j5old.sim.bly,
										ifelse(	presreg.court.outcome.bly.pivot=="smaller.shift" & above.bly==1, apply(  matrix(c(pres.sim.bly, (2*sen.pivot.sim.bly-j5old.sim.bly)),ncol=2) , 1, min),						
										ifelse(	presreg.court.outcome.bly.pivot=="smaller.shift" & above.bly==0, apply(  matrix(c(pres.sim.bly, (2*sen.pivot.sim.bly-j5old.sim.bly)),ncol=2) , 1, max),
										ifelse(	presreg.court.outcome.bly.pivot=="maxshift",pres.sim.bly, NA)))))			
									
# now define whether nominee is flip or shift within "smaller.shift"
# shift is if president is interior" to 2S-j -- can get his ideal point  in nearly court.outcome and mixed game
# flip  is if president is exterior to to 2S-j -- can only get S's reflection point  in nearly court.outcome and mixed game 3
#there are only about 300 simulations where shift is met						
	presreg.court.outcome.smaller.shift.type.bly.pivot <- 
			ifelse(above.bly==1 & 	presreg.court.outcome.bly.pivot=="smaller.shift"& pres.sim.bly <= 2*sen.pivot.sim.bly-j5old.sim.bly, "shift",
			ifelse(above.bly==1 & 	presreg.court.outcome.bly.pivot=="smaller.shift"& pres.sim.bly > 2*sen.pivot.sim.bly-j5old.sim.bly, "flip",
			ifelse(above.bly==0 & 	presreg.court.outcome.bly.pivot=="smaller.shift"& pres.sim.bly >= 2*sen.pivot.sim.bly-j5old.sim.bly, "shift",
			ifelse(above.bly==0 & 	presreg.court.outcome.bly.pivot=="smaller.shift"& pres.sim.bly < 2*sen.pivot.sim.bly-j5old.sim.bly, "flip", NA))))	
#NOW DO pos.taking REGIMES
		presreg.pos.taking.bly.pivot <- ifelse( above.bly==1 & (sen.regime.bly.pivot == "A" | sen.regime.bly.pivot=="B"), "gridlock",
									  	ifelse( above.bly==1 & (sen.regime.bly.pivot == "C" | sen.regime.bly.pivot=="D"), "smaller.shift",
										  ifelse( above.bly==0 & (sen.regime.bly.pivot == "C" | sen.regime.bly.pivot=="D"), "gridlock",
										  ifelse( above.bly==0 & (sen.regime.bly.pivot == "A" | sen.regime.bly.pivot=="B"), "smaller.shift", NA ))))
 #get point prediction for pos.taking
   pos.taking.point.pred.bly.pivot <- ifelse(presreg.pos.taking.bly.pivot=="gridlock", j5old.sim.bly,
										ifelse(	presreg.pos.taking.bly.pivot=="smaller.shift" & above.bly==1, apply(  matrix(c(pres.sim.bly, (2*sen.pivot.sim.bly-j5old.sim.bly)),ncol=2) , 1, min),						
										ifelse(	presreg.pos.taking.bly.pivot=="smaller.shift" & above.bly==0, apply(  matrix(c(pres.sim.bly, (2*sen.pivot.sim.bly-j5old.sim.bly)),ncol=2) , 1, max), NA)))																 
 	 presreg.pos.taking.smaller.shift.type.bly.pivot <- 
			ifelse(above.bly==1 & 	presreg.pos.taking.bly.pivot=="smaller.shift" & pres.sim.bly <= 2*sen.pivot.sim.bly-j5old.sim.bly, "shift",
			ifelse(above.bly==1 & 	presreg.pos.taking.bly.pivot=="smaller.shift" & pres.sim.bly > 2*sen.pivot.sim.bly-j5old.sim.bly, "flip",
			ifelse(above.bly==0 & 	presreg.pos.taking.bly.pivot=="smaller.shift" & pres.sim.bly >= 2*sen.pivot.sim.bly-j5old.sim.bly, "shift",
			ifelse(above.bly==0 & 	presreg.pos.taking.bly.pivot=="smaller.shift" & pres.sim.bly < 2*sen.pivot.sim.bly-j5old.sim.bly, "flip", NA))))	
	})

tapply(bailey.nomd.pivot$sen.regime.bly, bailey.nomd.pivot$nominee, table)
tapply(c(bailey.nomd.pivot$presreg.court.outcome.bly,bailey.nomd.pivot$sen.regime.bly), c(bailey.nomd.pivot$nominee,bailey.nomd.pivot$nominee) , table)
tapply(c(bailey.nomd.pivot$presreg.pos.taking.bly,bailey.nomd.pivot$sen.regime.bly), c(bailey.nomd.pivot$nominee,bailey.nomd.pivot$nominee) , table)

#make matrix with frequencies of all categories by nominee
library(zoo)
court.outcome.regime.dist.bly.pivot  <- data.frame(t(do.call(merge, lapply(tapply(bailey.nomd.pivot$presreg.court.outcome.bly.pivot, bailey.nomd.pivot$nominee, table), function(x) zoo(c(x), rownames(x))))))
court.outcome.regime.dist.bly.pivot <- within(court.outcome.regime.dist.bly.pivot, {
	nominee <- rownames(court.outcome.regime.dist.bly.pivot)
	court.outcome.regime.modal.bly.pivot  <- apply(court.outcome.regime.dist.bly.pivot ,1,which.max)
	court.outcome.regime.modal.bly.pivot <- ifelse(court.outcome.regime.modal.bly.pivot==1, "flip", ifelse(court.outcome.regime.modal.bly.pivot==2, "gridlock", 
											 ifelse(court.outcome.regime.modal.bly.pivot=="3", "maxshift", ifelse(court.outcome.regime.modal.bly.pivot==4,"restoring", NA))))
  court.outcome.regime.max.percent.bly.pivot <- apply(cbind(smaller.shift, gridlock, maxshift, restoring),1,max, na.rm=T)  / n.sims
})
head(court.outcome.regime.dist.bly.pivot)
hist(court.outcome.regime.dist.bly.pivot$court.outcome.regime.max.percent.bly.pivot)
table(court.outcome.regime.dist.bly.pivot$court.outcome.regime.max.percent.bly.pivot >.9)

#can we test flip or shift prediction?
table(court.outcome.regime.dist.bly.pivot$smaller.shift/n.sims>.5)#NO-

pos.taking.regime.dist.bly.pivot  <- data.frame(t(do.call(merge, lapply(tapply(bailey.nomd.pivot$presreg.pos.taking.bly.pivot, bailey.nomd.pivot$nominee, table), function(x) zoo(c(x), rownames(x))))))
pos.taking.regime.dist.bly.pivot <- within(pos.taking.regime.dist.bly.pivot, {
  nominee <- rownames(pos.taking.regime.dist.bly.pivot )
	pos.taking.regime.modal.bly.pivot  <- apply(pos.taking.regime.dist.bly.pivot ,1,which.max)
	pos.taking.regime.modal.bly.pivot <- ifelse(pos.taking.regime.modal.bly.pivot==1, "smaller.shift", ifelse(pos.taking.regime.modal.bly.pivot ==2, "gridlock",NA))
  pos.taking.regime.max.percent.bly.pivot <- apply(cbind(smaller.shift, gridlock),1,max, na.rm=T)  / n.sims
})
table(pos.taking.regime.dist.bly.pivot$pos.taking.regime.max.percent.bly.pivot>.8)

#separately, define "lower left" in terms of pres theory figure paper for testing 3rd robust prediction
#need to account for president
bailey.nomd.pivot <- within(bailey.nomd.pivot, {
	lower.left.bly.pivot <- ifelse( (above.bly ==1  & (sen.regime.bly.pivot == "A" | sen.regime.bly.pivot == "B") &  j.exit.sim.bly <= j5old.sim.bly )  | 
		(above.bly ==0  & (sen.regime.bly.pivot == "C" | sen.regime.bly.pivot == "D") &  j.exit.sim.bly > j5old.sim.bly) ,1,0)
})
prob.lower.left.bly.pivot <-   tapply (bailey.nomd.pivot$lower.left.bly, bailey.nomd.pivot$nominee,mean, na.omit=F)				

temp.df.pivot <- data.frame(nominee=court.outcome.regime.dist.bly.pivot$nominee,court.outcome.regime.modal.bly.pivot= court.outcome.regime.dist.bly.pivot$court.outcome.regime.modal.bly.pivot,
		court.outcome.regime.max.percent.bly.pivot = court.outcome.regime.dist.bly.pivot$court.outcome.regime.max.percent.bly.pivot ,
	 pos.taking.regime.modal.bly.pivot= pos.taking.regime.dist.bly.pivot$pos.taking.regime.modal.bly.pivot, 
		pos.taking.regime.max.percent.bly.pivot= pos.taking.regime.dist.bly.pivot$pos.taking.regime.max.percent.bly.pivot ,
		prob.lower.left.bly.pivot=prob.lower.left.bly.pivot )
	
#join to nom-data stats
bailey.nomd.pivot <- join(bailey.nomd.pivot,temp.df.pivot )
table(bailey.nomd.pivot$court.outcome.regime.modal.bly.pivot,bailey.nomd.pivot$pos.taking.regime.modal.bly.pivot )
head(bailey.nomd.pivot)

#COMPARE PREDICTIONS USING SENATE MEDIAN AND FILIBUSTER
	#TAKE MEANS BY NOMINEE
dw.preds <- ddply(nomd.pivot, .(nominee), summarize,
	ncourt.outcome.median = mean(ncourt.outcome.point.pred.dw),
	ncourt.outcome.pivot = mean(ncourt.outcome.point.pred.dw.pivot),
	pos.taking.median = mean(pos.taking.point.pred.dw),
	pos.taking.pivot = mean(pos.taking.point.pred.dw.pivot)
	)

bly.preds <- ddply(bailey.nomd.pivot, .(nominee), summarize,
	ncourt.outcome.median = mean(ncourt.outcome.point.pred.bly),
	ncourt.outcome.pivot = mean(ncourt.outcome.point.pred.bly.pivot),
	pos.taking.median = mean(pos.taking.point.pred.bly),
	pos.taking.pivot = mean(pos.taking.point.pred.bly.pivot)
	)

pdf("Figure_A4.pdf", height=8, width=8)

layout(rbind( c(9,5,6),c(7,1,3),c(8,2,4)), heights=c(.1,1,1), width=c(.1,1,1))
par(mar=c(4,4,1,2))
text.size <- .7
corr.label <- 1.2
label.text <- 1
axis.text <- 1.2
#NEARLY court.outcome MODEL
#DW
plot(dw.preds$ncourt.outcome.median, dw.preds$ncourt.outcome.pivot, xlim=c(-.5,.5), ylim=c(-.5,.5), pch=19, axes=F, xlab="", ylab="")
abline(0,1,lty=2)
axis(1, mgp=c(2,.5,.0), cex.axis=axis.text)
axis(2, mgp=c(2,.5,.0), las=1,cex.axis=axis.text)
mtext("Using median",1,line=2, cex=label.text)
mtext("Using pivot",2,line=3, cex=label.text,srt=90)
box()
#text(dw.preds$ncourt.outcome.median+.02, dw.preds$ncourt.outcome.pivot,dw.preds$nominee,adj=0, cex=.5)
text(.2,-.4, paste("p=", round(cor(dw.preds$ncourt.outcome.median, dw.preds$ncourt.outcome.pivot),2)), cex=corr.label)
#pos.taking MODEL
#DW
plot(dw.preds$pos.taking.median, dw.preds$pos.taking.pivot,xlim=c(-.5,.5), ylim=c(-.5,.5),  pch=19, axes=F, xlab="", ylab="")
abline(0,1,lty=2)
axis(1, mgp=c(2,.5,.0), cex.axis=axis.text)
axis(2, mgp=c(2,.5,.0), las=1,cex.axis=axis.text)
mtext("Using median",1,line=2, cex=label.text)
mtext("Using pivot",2,line=3, cex=label.text,srt=90)
box()
#text(dw.preds$pos.taking.median+.02, dw.preds$pos.taking.pivot,dw.preds$nominee,adj=0, cex=.5)
text(.2,-.4, paste("p=", round(cor(dw.preds$pos.taking.median, dw.preds$pos.taking.pivot),2)), cex=corr.label)

#BAILEY
plot(bly.preds$ncourt.outcome.median, bly.preds$ncourt.outcome.pivot, xlim=c(-1.2,1.2), ylim=c(-1.2,1.2), pch=19, axes=F, xlab="", ylab="")
#abline(lm(bly.preds$ncourt.outcome.pivot~bly.preds$ncourt.outcome.median))
abline(0,1,lty=2)
axis(1, mgp=c(2,.5,.0), cex.axis=axis.text)
axis(2, mgp=c(2,.5,.0), las=1,cex.axis=axis.text)
mtext("Using median",1,line=2, cex=label.text)
mtext("Using pivot",2,line=3, cex=label.text,srt=90)
box()
text(.5,-.9, paste("p=", round(cor(bly.preds$ncourt.outcome.median, bly.preds$ncourt.outcome.pivot),2)), cex=corr.label)

plot(bly.preds$pos.taking.median, bly.preds$pos.taking.pivot, xlim=c(-1.2,1.2), ylim=c(-1.2,1.2), pch=19, axes=F, xlab="", ylab="")
#abline(lm(bly.preds$pos.taking.pivot~bly.preds$pos.taking.median))
abline(0,1,lty=2)
axis(1, mgp=c(2,.5,.0), cex.axis=axis.text)
axis(2, mgp=c(2,.5,.0), las=1,cex.axis=axis.text)
mtext("Using median",1,line=2, cex=label.text)
mtext("Using pivot",2,line=3, cex=label.text,srt=90)
box()
text(.5,-.9, paste("p=", round(cor(bly.preds$pos.taking.median, bly.preds$pos.taking.pivot),2)), cex=corr.label)

#text(bly.preds$pos.taking.median+.02, bly.preds$pos.taking.pivot,bly.preds$nominee,adj=0, cex=.5)
par(mar=c(0,0,0,0))
plot(0,0,type="n", xlab="", ylab="",main="", axes=F)
text(0,-.5, "NOMINATE", xpd=NA, font=2, cex=1.5)
plot(0,0,type="n", xlab="", ylab="",main="", axes=F)
text(0,-.5, "BAILEY", xpd=NA, font=2, cex=1.5)
plot(0,0,type="n", xlab="", ylab="",main="", axes=F)
text(0,.1, "Nearly court-outcome based", xpd=NA, font=2, cex=2, srt=90)
plot(0,0,type="n", xlab="", ylab="",main="", axes=F)
text(0,.1, "Position-taking senators", xpd=NA, font=2, cex=2, srt=90)

dev.off()

##########
#TEST MEDIAN-locked prediction
##########
#N: take nominations where predicted to be in gridlock across all regimes
#compare actual versus j5ol
#NOMINATE
nomd.pivot$gridlock.robust.pivot <- ifelse(nomd.pivot$presreg.court.outcome.dw.pivot=="gridlock" & nomd.pivot$presreg.pos.taking.dw.pivot=="gridlock","gridlock","not gridlock")
gridlock.robust.dw.pivot  <- data.frame(t(do.call(merge, lapply(tapply(nomd.pivot$gridlock.robust , nomd.pivot$nominee, table), function(x) zoo(c(x), rownames(x))))))
gridlock.robust.dw.pivot$gridlock <- ifelse(is.na(gridlock.robust.dw.pivot$gridlock), 0, gridlock.robust.dw.pivot$gridlock)
gridlock.robust.dw.pivot$not.gridlock <- ifelse(is.na(gridlock.robust.dw.pivot$not.gridlock), 0, gridlock.robust.dw.pivot$not.gridlock)

gridlock.robust.dw.pivot <- within(gridlock.robust.dw.pivot, {
	nominee <- rownames(gridlock.robust.dw.pivot)
	gridlock.robust.modal.dw.pivot  <- apply(gridlock.robust.dw.pivot ,1,which.max)
	gridlock.robust.modal.dw.pivot <- ifelse(gridlock.robust.modal.dw.pivot==2, "not gridlock", ifelse(gridlock.robust.modal.dw.pivot==1, "gridlock", NA))
  gridlock.robust.max.percent.dw.pivot <- apply(cbind(gridlock),1,max, na.rm=T)  / n.sims
})
head(gridlock.robust.dw.pivot)
hist(gridlock.robust.dw.pivot$gridlock.robust.max.percent.dw.pivot)
table(gridlock.robust.dw.pivot$gridlock.robust.max.percent.dw.pivot>.5)
nomd.pivot <- join(nomd.pivot,gridlock.robust.dw.pivot)

#Bailey
bailey.nomd.pivot$gridlock.robust.pivot <- ifelse(bailey.nomd.pivot$presreg.court.outcome.bly.pivot=="gridlock" & bailey.nomd.pivot$presreg.pos.taking.bly.pivot=="gridlock","gridlock","not gridlock")
gridlock.robust.bly.pivot  <- data.frame(t(do.call(merge, lapply(tapply(bailey.nomd.pivot$gridlock.robust , bailey.nomd.pivot$nominee, table), function(x) zoo(c(x), rownames(x))))))
gridlock.robust.bly.pivot$gridlock <- ifelse(is.na(gridlock.robust.bly.pivot$gridlock), 0, gridlock.robust.bly.pivot$gridlock)
gridlock.robust.bly.pivot$not.gridlock <- ifelse(is.na(gridlock.robust.bly.pivot$not.gridlock), 0, gridlock.robust.bly.pivot$not.gridlock)

gridlock.robust.bly.pivot <- within(gridlock.robust.bly.pivot, {
	nominee <- rownames(gridlock.robust.bly.pivot)
	gridlock.robust.modal.bly.pivot  <- apply(gridlock.robust.bly.pivot ,1,which.max)
	gridlock.robust.modal.bly.pivot <- ifelse(gridlock.robust.modal.bly.pivot==2, "not gridlock", ifelse(gridlock.robust.modal.bly.pivot==1, "gridlock", NA))
  gridlock.robust.max.percent.bly.pivot <- apply(cbind(gridlock),1,max, na.rm=T)  / n.sims
})
head(gridlock.robust.bly.pivot)
hist(gridlock.robust.bly.pivot$gridlock.robust.max.percent.bly.pivot)
table(gridlock.robust.bly.pivot$gridlock.robust.max.percent.bly.pivot>.5)
bailey.nomd.pivot <- join(bailey.nomd.pivot,gridlock.robust.bly.pivot)

#regress: completely constrained 
threshold.hi <-.5

#MAKE GRAPH

pdf("Figure_A6.pdf", height=10, width=12)

x.axis <-1.2
y.axis.text <-1.1
label.text <- 1.2

#dw.pivot
ok <- nomd.pivot$gridlock.robust.max.percent.dw.pivot >= threshold.hi
temp.grid.data.dw.pivot <- subset(nomd.pivot, ok)#make temp data
#calculate confidence intervals and median estimate

temp.grid.data.dw.pivot$nom.v.old <- temp.grid.data.dw.pivot$nom.proj.sim.dw-temp.grid.data.dw.pivot$j5old.sim.dw

nom.v.old.lo.dw.pivot <- tapply(temp.grid.data.dw.pivot$nom.v.old, temp.grid.data.dw.pivot$nominee, quantile,probs=c(lo.ci))
nom.v.old.med.dw.pivot <-  tapply(temp.grid.data.dw.pivot$nom.v.old, temp.grid.data.dw.pivot$nominee, quantile,probs=c(.5))
nom.v.old.hi.dw.pivot <- tapply(temp.grid.data.dw.pivot$nom.v.old, temp.grid.data.dw.pivot$nominee, quantile,probs=c(hi.ci))
nom.v.old.mean.prob.dw.pivot <- tapply(temp.grid.data.dw.pivot$nom.v.old > 0, temp.grid.data.dw.pivot$nominee,mean)

dem.pres.temp <- ifelse(tapply(temp.grid.data.dw.pivot$pres.sim.dw<0,temp.grid.data.dw.pivot$nominee, mean)>0,1,0)
confirmed.temp <-  tapply(temp.grid.data.dw.pivot$confirmed > 0, temp.grid.data.dw.pivot$nominee,mean)
temp.df <- data.frame(nominee=unique(temp.grid.data.dw.pivot$nominee.title.dw),
	nom.v.old.lo.dw.pivot,nom.v.old.med.dw.pivot,nom.v.old.hi.dw.pivot, nom.v.old.mean.prob.dw.pivot, dem.pres.temp, confirmed.temp)
temp.df <- temp.df[order(temp.df$dem.pres.temp,nom.v.old.med.dw.pivot ),]
	
y <- 1:nrow(temp.df)
x <-temp.df$nom.v.old.med.dw.pivot
x.lo <- temp.df$nom.v.old.lo.dw.pivot
x.hi <- temp.df$nom.v.old.hi.dw.pivot
nominee.sorted <- temp.df$nominee
confirmed.sorted <- temp.df$confirmed.temp
prob.text <- .8

par(mfrow=c(1,2),mar=c(3,6.5,.5,1))
plot(x,y , xlim=c(-.5,.5), pch=19, col="black", xlab="", ylab="",axes=F, type="n")
polygon(x=c(par()$usr[1],par()$usr[1], par()$usr[2],par()$usr[2]), y=c(11.5,par()$usr[4]-.5,par()$usr[4]-.5,11.5), col="gray90")
segments(x.lo, y, x.hi, y)
points(x[confirmed.sorted==1], y[confirmed.sorted==1], pch=19)
points(x[confirmed.sorted==0], y[confirmed.sorted==0], pch=21, bg="white")

axis(1, at = seq(-1,1.4,.2), 
     labels=c("-1","-.8","-.6",  "-.4", "-.2", 0, ".2", ".4", ".6",".8","1","1.2","1.4"), mgp=c(2,.5,0), cex.axis=x.axis)
axis(2, at =y , label=nominee.sorted, mgp=c(2,.5,0), cex.axis=y.axis.text, las=2)
segments(0,par()$usr[4]-.5,0,par()$usr[3], col="gray80")

mtext("Nominee-old median", 1, line=2,cex=label.text)
mtext("NOMINATE",3, line=-1, cex=label.text, font=2)

text(-.4, 18, "Democratic\nappointees",cex=label.text, font=3, srt=90)

#Bailey
ok <- bailey.nomd.pivot$gridlock.robust.max.percent.bly >= threshold.hi
temp.grid.data.bly.pivot <- subset(bailey.nomd.pivot, ok)#make temp data
#calculate confidence intervals and median estimate
temp.grid.data.bly.pivot$nom.v.old <- temp.grid.data.bly.pivot$nom.proj.sim.bly-temp.grid.data.bly.pivot$j5old.sim.bly 

nom.v.old.lo.bly <- tapply(temp.grid.data.bly.pivot$nom.v.old, temp.grid.data.bly.pivot$nominee, quantile,probs=c(lo.ci))
nom.v.old.med.bly <-  tapply(temp.grid.data.bly.pivot$nom.v.old, temp.grid.data.bly.pivot$nominee, quantile,probs=c(.5))
nom.v.old.hi.bly <- tapply(temp.grid.data.bly.pivot$nom.v.old, temp.grid.data.bly.pivot$nominee, quantile,probs=c(hi.ci))
nom.v.old.mean.prob.bly <- tapply(temp.grid.data.bly.pivot$nom.v.old > 0, temp.grid.data.bly.pivot$nominee,mean)

dem.pres.temp <- ifelse( tapply(temp.grid.data.bly.pivot$pres.sim.bly<0,temp.grid.data.bly.pivot$nominee, mean) > .5,1,0)#need cutoff of .5 to separate
confirmed.temp <-  tapply(temp.grid.data.bly.pivot$confirmed > 0, temp.grid.data.bly.pivot$nominee,mean)
temp.df <- data.frame(nominee=unique(temp.grid.data.bly.pivot$nominee.title.bly),
	nom.v.old.lo.bly,nom.v.old.med.bly,nom.v.old.hi.bly, nom.v.old.mean.prob.bly, dem.pres.temp, confirmed.temp)
temp.df <- temp.df[order(temp.df$dem.pres.temp,nom.v.old.med.bly ),]
	
y <- 1:nrow(temp.df)
x <-temp.df$nom.v.old.med.bly
x.lo <- temp.df$nom.v.old.lo.bly
x.hi <- temp.df$nom.v.old.hi.bly
nominee.sorted <- temp.df$nominee
confirmed.sorted <- temp.df$confirmed.temp
prob.text <- .8

#confirmed.sorted <- bailey.nom.level.stats.pivot$confirmed[order(bailey.nom.level.stats.pivot$agg.mistake.court.outcome.distance.med.bly)]
plot(x,y , xlim=c(-1.4,1.4), pch=19, col="black", xlab="", ylab="",axes=F, type="n")
polygon(x=c(par()$usr[1],par()$usr[1], par()$usr[2],par()$usr[2]), y=c(11.5,par()$usr[4]-.4,par()$usr[4]-.4,11.5), col="gray90")
segments(x.lo, y, x.hi, y)
points(x[confirmed.sorted==1], y[confirmed.sorted==1], pch=19)
points(x[confirmed.sorted==0], y[confirmed.sorted==0], pch=21, bg="white")
axis(1, at = seq(-1.4,1.4,.2), 
     labels=c("-1.4","-1.2","-1","-.8","-.6",  "-.4", "-.2", 0, ".2", ".4", ".6",".8","1","1.2","1.4"), mgp=c(2,.5,0), cex.axis=x.axis)
axis(2, at =y , label=nominee.sorted, mgp=c(2,.5,0), cex.axis=y.axis.text, las=2)
segments(0,par()$usr[4]-.5,0,par()$usr[3], col="gray80")
mtext("Nominee-old median", 1, line=2,cex=label.text)
mtext("Bailey",3, line=-1, cex=label.text, font=2)
text(-1.2, 13, "Democratic\nappointees",cex=label.text, font=3, srt=90)

dev.off()
#how many nominees include zero?
table(nom.v.old.lo.dw<0 & nom.v.old.hi.dw>0)
table(nom.v.old.lo.bly<0 & nom.v.old.hi.bly>0)

#PRESIDENTIAL REGRESSIONS USING PIVOT (Table A-3)
########################
#DW
########################
nomd.pivot <- within (nomd.pivot, {
	#court.outcome MODELS
		#CREATE DUMMIES
	court.outcome.gridlock.dw.pivot <- ifelse(presreg.court.outcome.dw.pivot=="gridlock", 1,0)
	court.outcome.flip.dw.pivot <- ifelse(presreg.court.outcome.smaller.shift.type.dw.pivot=="flip", 1,0)
	court.outcome.flip.dw.pivot <- ifelse( is.na(court.outcome.flip.dw.pivot), 0,court.outcome.flip.dw.pivot)
	court.outcome.shift.dw.pivot <- ifelse(presreg.court.outcome.smaller.shift.type.dw.pivot=="shift", 1,0)
	court.outcome.shift.dw.pivot <- ifelse( is.na(court.outcome.shift.dw.pivot), 0,court.outcome.shift.dw.pivot)
	court.outcome.restoring.dw.pivot <- ifelse(presreg.court.outcome.dw.pivot=="restoring", 1,0)
	court.outcome.maxshift.dw.pivot <- ifelse(presreg.court.outcome.dw.pivot=="maxshift", 1,0)
	court.outcome.p.combined.dw.pivot <-ifelse(	court.outcome.restoring.dw.pivot==1 | 	court.outcome.maxshift.dw.pivot==1 |	court.outcome.shift.dw.pivot== 1 , 1,0)
	court.outcome.p.combined.dw.pivot <- ifelse( is.na(court.outcome.p.combined.dw.pivot), 0,court.outcome.p.combined.dw.pivot)
	#for combining regions
	court.outcome.not.p.combined.dw.pivot <- ifelse(	court.outcome.p.combined.dw.pivot==0,1,0)
	court.outcome.not.gridlock.dw.pivot <- ifelse(	court.outcome.gridlock.dw.pivot==0,1,0)
	court.outcome.not.flip.dw.pivot <- ifelse(court.outcome.flip.dw.pivot==0,1,0)
	#pos.taking MODELS
	pos.taking.gridlock.dw.pivot <- ifelse(presreg.pos.taking.dw.pivot=="gridlock", 1,0)
	pos.taking.flip.dw.pivot     <- ifelse(presreg.pos.taking.smaller.shift.type.dw.pivot== "flip", 1,0)
	pos.taking.flip.dw.pivot     <- ifelse( is.na(pos.taking.flip.dw.pivot ), 0,pos.taking.flip.dw.pivot )
	pos.taking.shift.dw.pivot    <- ifelse(presreg.pos.taking.smaller.shift.type.dw.pivot== "shift", 1,0)
	pos.taking.shift.dw.pivot    <- ifelse( is.na(pos.taking.shift.dw.pivot ), 0,pos.taking.shift.dw.pivot )
	pos.taking.not.gridlock.dw.pivot <- ifelse(	pos.taking.gridlock.dw.pivot==0,1,0)
	pos.taking.not.flip.dw.pivot    <- ifelse(	pos.taking.flip.dw.pivot==0,1,0)
	pos.taking.not.shift.dw.pivot   <- ifelse(	pos.taking.shift.dw.pivot==0,1,0)
	#Senate indifferen      ce POINT: this applies across models
	sen.ind.point.dw.pivot   <- 2*sen.pivot.sim.dw-j5old.sim.dw	 
})

#confirm implicated regimes are different across median/FP
plot(nomd.pivot$sen.ind.point.dw.pivot, nomd$sen.ind.point.dw)
table(nomd.pivot$	court.outcome.gridlock.dw.pivot, nomd$	court.outcome.gridlock.dw)
table(nomd.pivot$	court.outcome.flip.dw.pivot, nomd$		court.outcome.flip.dw)
table(nomd.pivot$		pos.taking.gridlock.dw.pivot , nomd$	pos.taking.gridlock.dw)
table(nomd.pivot$	court.outcome.p.combined.dw.pivot, nomd$		court.outcome.p.combined.dw)


lo <-.025
hi <- .975
#ncourt.outcome model#
#Model 1: just theory predictions
court.outcome.model.pivot.1.dw.intercept <- rep(NA, n.sims)
court.outcome.model.pivot.1.dw.gridlock <- rep(NA, n.sims)
court.outcome.model.pivot.1.dw.p.predicted <- rep(NA, n.sims)
court.outcome.model.pivot.1.dw.flip <- rep(NA, n.sims)
court.outcome.model.pivot.1.dw.rsquared <- rep(NA, n.sims)

for (i in 1:n.sims){
	ok <- nomd.pivot$simulation==i
  foo <-	lm(nom.proj.sim.dw ~  	court.outcome.gridlock.dw.pivot:j5old.sim.dw + court.outcome.p.combined.dw.pivot:pres.sim.dw +	court.outcome.flip.dw.pivot:sen.ind.point.dw ,	
		data=nomd.pivot,subset=ok)
	sim.foo <- sim(foo,n.sims=1)
	court.outcome.model.pivot.1.dw.intercept[i] <- coef(sim.foo)[1]
	court.outcome.model.pivot.1.dw.gridlock[i] <- coef(sim.foo)[2]
	court.outcome.model.pivot.1.dw.p.predicted[i] <- coef(sim.foo)[3]
	court.outcome.model.pivot.1.dw.flip[i] <- coef(sim.foo)[4]
	court.outcome.model.pivot.1.dw.rsquared[i] <- summary(foo)$r.squared
}		

#MODEL 2: now in non-theory predictions
court.outcome.model.pivot.2.dw.intercept <- rep(NA, n.sims)
court.outcome.model.pivot.2.dw.not.gridlock <- rep(NA, n.sims)
court.outcome.model.pivot.2.dw.not.p.predicted <- rep(NA, n.sims)
court.outcome.model.pivot.2.dw.not.flip <- rep(NA, n.sims)
court.outcome.model.pivot.2.dw.rsquared <- rep(NA, n.sims)

for (i in 1:n.sims){
	ok <- nomd.pivot$simulation==i
  foo <-	lm(nom.proj.sim.dw ~  	
		+ 	court.outcome.not.gridlock.dw.pivot:j5old.sim.dw + court.outcome.not.p.combined.dw.pivot:pres.sim.dw +	court.outcome.not.flip.dw.pivot:sen.ind.point.dw
		, data=nomd.pivot,subset=ok)
	sim.foo <- sim(foo,n.sims=1)
	court.outcome.model.pivot.2.dw.intercept[i]    <- coef(sim.foo)[1]
	court.outcome.model.pivot.2.dw.not.gridlock[i]    <- coef(sim.foo)[2]
	court.outcome.model.pivot.2.dw.not.p.predicted[i] <- coef(sim.foo)[3]
	court.outcome.model.pivot.2.dw.not.flip[i]        <- coef(sim.foo)[4]	
	court.outcome.model.pivot.2.dw.rsquared[i] <- summary(foo)$r.squared
}		
############################################
######################
#pos.taking model
######################
######################
#1) just theory predictions
pos.taking.model.pivot.1.dw.intercept <- rep(NA, n.sims)
pos.taking.model.pivot.1.dw.gridlock <- rep(NA, n.sims)
pos.taking.model.pivot.1.dw.shift <- rep(NA, n.sims)
pos.taking.model.pivot.1.dw.flip <- rep(NA, n.sims)
pos.taking.model.pivot.1.dw.rsquared <- rep(NA, n.sims)
	
for (i in 1:n.sims){
	ok <- nomd.pivot$simulation==i
	#pos.taking
 foo <-	lm(nom.proj.sim.dw ~  	pos.taking.gridlock.dw.pivot:j5old.sim.dw + pos.taking.shift.dw.pivot:pres.sim.dw + pos.taking.flip.dw.pivot:sen.ind.point.dw,
		 	data=nomd.pivot,subset=ok)
	#display(foo)
	sim.foo <- sim(foo,n.sims=1)
	pos.taking.model.pivot.1.dw.intercept[i] <- coef(sim.foo)[1]
	pos.taking.model.pivot.1.dw.gridlock[i] <- coef(sim.foo)[2]
	pos.taking.model.pivot.1.dw.shift[i] <- coef(sim.foo)[3]
	pos.taking.model.pivot.1.dw.flip[i] <- coef(sim.foo)[4]
	pos.taking.model.pivot.1.dw.rsquared[i] <- summary(foo)$r.squared
}		

#2) now non-theory specific predictgions
pos.taking.model.pivot.2.dw.intercept <- rep(NA, n.sims)
pos.taking.model.pivot.2.dw.not.gridlock <- rep(NA, n.sims)
pos.taking.model.pivot.2.dw.not.shift <- rep(NA, n.sims)
pos.taking.model.pivot.2.dw.not.flip <- rep(NA, n.sims)
pos.taking.model.pivot.2.dw.rsquared <- rep(NA, n.sims)

for (i in 1:n.sims){
	ok <- nomd.pivot$simulation==i
	#pos.taking
 foo <-	lm(nom.proj.sim.dw ~  	pos.taking.not.gridlock.dw.pivot:j5old.sim.dw + pos.taking.not.shift.dw.pivot:pres.sim.dw + pos.taking.not.flip.dw.pivot:sen.ind.point.dw, data=nomd.pivot,subset=ok)
	sim.foo <- sim(foo,n.sims=1)
	pos.taking.model.pivot.2.dw.intercept[i] <- coef(foo)[1]
	pos.taking.model.pivot.2.dw.not.gridlock[i] <-  coef(foo)[2]
	pos.taking.model.pivot.2.dw.not.shift[i] <- coef(foo)[3]
  pos.taking.model.pivot.2.dw.not.flip[i]  <-	 coef(foo)[4]
	pos.taking.model.pivot.2.dw.rsquared[i] <- summary(foo)$r.squared
}

#####################################
#BAILEY
#####################################
bailey.nomd.pivot <- within (bailey.nomd.pivot, {
	#court.outcome MODELS
		#CREATE DUMMIES
	court.outcome.gridlock.bly.pivot <- ifelse(presreg.court.outcome.bly.pivot=="gridlock", 1,0)
	court.outcome.flip.bly.pivot <- ifelse(presreg.court.outcome.smaller.shift.type.bly.pivot=="flip", 1,0)
	court.outcome.flip.bly.pivot <- ifelse( is.na(court.outcome.flip.bly.pivot), 0,court.outcome.flip.bly.pivot)
	court.outcome.shift.bly.pivot <- ifelse(presreg.court.outcome.smaller.shift.type.bly.pivot=="shift", 1,0)
	court.outcome.shift.bly.pivot <- ifelse( is.na(court.outcome.shift.bly.pivot), 0,court.outcome.shift.bly.pivot)
	court.outcome.restoring.bly.pivot <- ifelse(presreg.court.outcome.bly.pivot=="restoring", 1,0)
	court.outcome.maxshift.bly.pivot <- ifelse(presreg.court.outcome.bly.pivot=="maxshift", 1,0)
	court.outcome.p.combined.bly.pivot <-ifelse(	court.outcome.restoring.bly.pivot==1 | 	court.outcome.maxshift.bly.pivot==1 |	court.outcome.shift.bly.pivot== 1 , 1,0)
	court.outcome.p.combined.bly.pivot <- ifelse( is.na(court.outcome.p.combined.bly.pivot), 0,court.outcome.p.combined.bly.pivot)
	#for combining regions
	court.outcome.not.p.combined.bly.pivot <- ifelse(	court.outcome.p.combined.bly.pivot==0,1,0)
	court.outcome.not.gridlock.bly.pivot <- ifelse(	court.outcome.gridlock.bly.pivot==0,1,0)
	court.outcome.not.flip.bly.pivot <- ifelse(court.outcome.flip.bly.pivot==0,1,0)
	#pos.taking MODELS
	pos.taking.gridlock.bly.pivot <- ifelse(presreg.pos.taking.bly.pivot=="gridlock", 1,0)
	pos.taking.flip.bly.pivot     <- ifelse(presreg.pos.taking.smaller.shift.type.bly.pivot== "flip", 1,0)
	pos.taking.flip.bly.pivot     <- ifelse( is.na(pos.taking.flip.bly.pivot ), 0,pos.taking.flip.bly.pivot )
	pos.taking.shift.bly.pivot    <- ifelse(presreg.pos.taking.smaller.shift.type.bly.pivot== "shift", 1,0)
	pos.taking.shift.bly.pivot    <- ifelse( is.na(pos.taking.shift.bly.pivot ), 0,pos.taking.shift.bly.pivot )
	pos.taking.not.gridlock.bly.pivot <- ifelse(	pos.taking.gridlock.bly.pivot==0,1,0)
	pos.taking.not.flip.bly.pivot    <- ifelse(	pos.taking.flip.bly.pivot==0,1,0)
	pos.taking.not.shift.bly.pivot   <- ifelse(	pos.taking.shift.bly.pivot==0,1,0)
	#Senate indifferen      ce POINT: this applies across models
	sen.ind.point.bly.pivot   <- 2*sen.pivot.sim.bly-j5old.sim.bly	 
})

lo <-.025
hi <- .975
#ncourt.outcome model#
#Model 1: just theory predictions
court.outcome.model.pivot.1.bly.intercept <- rep(NA, n.sims)
court.outcome.model.pivot.1.bly.gridlock <- rep(NA, n.sims)
court.outcome.model.pivot.1.bly.p.predicted <- rep(NA, n.sims)
court.outcome.model.pivot.1.bly.flip <- rep(NA, n.sims)
court.outcome.model.pivot.1.bly.rsquared <- rep(NA, n.sims)

for (i in 1:n.sims){
	ok <- bailey.nomd.pivot$simulation==i
  foo <-	lm(nom.proj.sim.bly ~  	court.outcome.gridlock.bly.pivot:j5old.sim.bly + court.outcome.p.combined.bly.pivot:pres.sim.bly +	court.outcome.flip.bly.pivot:sen.ind.point.bly ,	
		data=bailey.nomd.pivot,subset=ok)
	sim.foo <- sim(foo,n.sims=1)
	court.outcome.model.pivot.1.bly.intercept[i] <- coef(sim.foo)[1]
	court.outcome.model.pivot.1.bly.gridlock[i] <- coef(sim.foo)[2]
	court.outcome.model.pivot.1.bly.p.predicted[i] <- coef(sim.foo)[3]
	court.outcome.model.pivot.1.bly.flip[i] <- coef(sim.foo)[4]
	court.outcome.model.pivot.1.bly.rsquared[i] <- summary(foo)$r.squared
}		

#MODEL 2: now in non-theory predictions
court.outcome.model.pivot.2.bly.intercept <- rep(NA, n.sims)
court.outcome.model.pivot.2.bly.not.gridlock <- rep(NA, n.sims)
court.outcome.model.pivot.2.bly.not.p.predicted <- rep(NA, n.sims)
court.outcome.model.pivot.2.bly.not.flip <- rep(NA, n.sims)
court.outcome.model.pivot.2.bly.rsquared <- rep(NA, n.sims)

for (i in 1:n.sims){
	ok <- bailey.nomd.pivot$simulation==i
  foo <-	lm(nom.proj.sim.bly ~  	
		+ 	court.outcome.not.gridlock.bly.pivot:j5old.sim.bly + court.outcome.not.p.combined.bly.pivot:pres.sim.bly +	court.outcome.not.flip.bly.pivot:sen.ind.point.bly
		, data=bailey.nomd.pivot,subset=ok)
	sim.foo <- sim(foo,n.sims=1)
	court.outcome.model.pivot.2.bly.intercept[i]    <- coef(sim.foo)[1]
	court.outcome.model.pivot.2.bly.not.gridlock[i]    <- coef(sim.foo)[2]
	court.outcome.model.pivot.2.bly.not.p.predicted[i] <- coef(sim.foo)[3]
	court.outcome.model.pivot.2.bly.not.flip[i]        <- coef(sim.foo)[4]	
	court.outcome.model.pivot.2.bly.rsquared[i] <- summary(foo)$r.squared
}		
############################################
######################
#pos.taking model
######################
######################
#1) just theory predictions
pos.taking.model.pivot.1.bly.intercept <- rep(NA, n.sims)
pos.taking.model.pivot.1.bly.gridlock <- rep(NA, n.sims)
pos.taking.model.pivot.1.bly.shift <- rep(NA, n.sims)
pos.taking.model.pivot.1.bly.flip <- rep(NA, n.sims)
pos.taking.model.pivot.1.bly.rsquared <- rep(NA, n.sims)
	
for (i in 1:n.sims){
	ok <- bailey.nomd.pivot$simulation==i
	#pos.taking
 foo <-	lm(nom.proj.sim.bly ~  	pos.taking.gridlock.bly.pivot:j5old.sim.bly + pos.taking.shift.bly.pivot:pres.sim.bly + pos.taking.flip.bly.pivot:sen.ind.point.bly,
		 	data=bailey.nomd.pivot,subset=ok)
	#display(foo)
	sim.foo <- sim(foo,n.sims=1)
	pos.taking.model.pivot.1.bly.intercept[i] <- coef(sim.foo)[1]
	pos.taking.model.pivot.1.bly.gridlock[i] <- coef(sim.foo)[2]
	pos.taking.model.pivot.1.bly.shift[i] <- coef(sim.foo)[3]
	pos.taking.model.pivot.1.bly.flip[i] <- coef(sim.foo)[4]
	pos.taking.model.pivot.1.bly.rsquared[i] <- summary(foo)$r.squared
}		

#2) now non-theory specific predictgions
pos.taking.model.pivot.2.bly.intercept <- rep(NA, n.sims)
pos.taking.model.pivot.2.bly.not.gridlock <- rep(NA, n.sims)
pos.taking.model.pivot.2.bly.not.shift <- rep(NA, n.sims)
pos.taking.model.pivot.2.bly.not.flip <- rep(NA, n.sims)
pos.taking.model.pivot.2.bly.rsquared <- rep(NA, n.sims)

for (i in 1:n.sims){
	ok <- bailey.nomd.pivot$simulation==i
	#pos.taking
 foo <-	lm(nom.proj.sim.bly ~  	pos.taking.not.gridlock.bly.pivot:j5old.sim.bly + pos.taking.not.shift.bly.pivot:pres.sim.bly + pos.taking.not.flip.bly.pivot:sen.ind.point.bly, data=bailey.nomd.pivot,subset=ok)
	sim.foo <- sim(foo,n.sims=1)
	pos.taking.model.pivot.2.bly.intercept[i] <- coef(foo)[1]
	pos.taking.model.pivot.2.bly.not.gridlock[i] <-  coef(foo)[2]
	pos.taking.model.pivot.2.bly.not.shift[i] <- coef(foo)[3]
  pos.taking.model.pivot.2.bly.not.flip[i]  <-	 coef(foo)[4]
	pos.taking.model.pivot.2.bly.rsquared[i] <- summary(foo)$r.squared
}

#REPORT RESULTS FOR TABLE
#NEARLY court.outcome
	#CORRECT MODELS
#DW
quantile(court.outcome.model.pivot.1.dw.intercept , c(lo,.5,hi))
quantile(court.outcome.model.pivot.1.dw.gridlock, c(lo,.5,hi))
quantile(	court.outcome.model.pivot.1.dw.p.predicted, c(lo,.5,hi))
quantile(	court.outcome.model.pivot.1.dw.flip , c(lo,.5,hi), na.rm=T)
mean(	court.outcome.model.pivot.1.dw.rsquared)
#BLY
quantile(court.outcome.model.pivot.1.bly.intercept , c(lo,.5,hi))
quantile(court.outcome.model.pivot.1.bly.gridlock, c(lo,.5,hi))
quantile(	court.outcome.model.pivot.1.bly.p.predicted, c(lo,.5,hi))
quantile(	court.outcome.model.pivot.1.bly.flip , c(lo,.5,hi), na.rm=T)
mean(	court.outcome.model.pivot.1.bly.rsquared)
#PLACEBO MODELS
#DW
quantile(court.outcome.model.pivot.2.dw.intercept , c(lo,.5,hi))
quantile(court.outcome.model.pivot.2.dw.not.gridlock, c(lo,.5,hi))
quantile(court.outcome.model.pivot.2.dw.not.p.predicted, c(lo,.5,hi))
quantile(court.outcome.model.pivot.2.dw.not.flip , c(lo,.5,hi), na.rm=T)
mean(	court.outcome.model.pivot.2.dw.rsquared)
#BLY
quantile(court.outcome.model.pivot.2.bly.intercept , c(lo,.5,hi))
quantile(court.outcome.model.pivot.2.bly.not.gridlock, c(lo,.5,hi))
quantile(court.outcome.model.pivot.2.bly.not.p.predicted, c(lo,.5,hi))
quantile(court.outcome.model.pivot.2.bly.not.flip , c(lo,.5,hi), na.rm=T)
mean(	court.outcome.model.pivot.2.bly.rsquared)

#pos.taking
#CORRECT
#DW
quantile(	pos.taking.model.pivot.1.dw.intercept, c(lo,.5,hi))
quantile(	pos.taking.model.pivot.1.dw.gridlock, c(lo,.5,hi))
quantile(	pos.taking.model.pivot.1.dw.shift, c(lo,.5,hi), na.rm=T)
quantile(	pos.taking.model.pivot.1.dw.flip, c(lo,.5,hi), na.rm=T)
mean(	pos.taking.model.pivot.1.dw.rsquared)
#Bailey
quantile(	pos.taking.model.pivot.1.bly.intercept, c(lo,.5,hi))
quantile(	pos.taking.model.pivot.1.bly.gridlock, c(lo,.5,hi))
quantile(	pos.taking.model.pivot.1.bly.shift, c(lo,.5,hi))
quantile(	pos.taking.model.pivot.1.bly.flip, c(lo,.5,hi), na.rm=T)
mean(	pos.taking.model.pivot.1.bly.rsquared)
#PLACEBO
quantile( pos.taking.model.pivot.2.dw.intercept, c(lo,.5,hi))
quantile(	pos.taking.model.pivot.2.dw.not.gridlock, c(lo,.5,hi))
quantile(	pos.taking.model.pivot.2.dw.not.shift, c(lo,.5,hi))
quantile(	pos.taking.model.pivot.2.dw.not.flip, c(lo,.5,hi), na.rm=T)
mean(	pos.taking.model.pivot.2.dw.rsquared)

quantile( pos.taking.model.pivot.2.bly.intercept, c(lo,.5,hi))
quantile(	pos.taking.model.pivot.2.bly.not.gridlock, c(lo,.5,hi))
quantile(	pos.taking.model.pivot.2.bly.not.shift, c(lo,.5,hi))
quantile(	pos.taking.model.pivot.2.bly.not.flip, c(lo,.5,hi), na.rm=T)
mean(	pos.taking.model.pivot.2.bly.rsquared)



